Maths for Intelligent Systems
Marc ToussaintApril, 2022

This script is primarily based on a lecture I gave 2015-2019 in Stuttgart. The current
form also integrates notes and exercises from the Optimization lecture, and a little
material from my Robotics and ML lectures. The first full version was from 2019, sincethen I occasionally update it.

Contents
1 Speaking Maths

5

1.1

Describing systems. . . . . . . . . . . . . . . . . . . . . . . . . . . .

5

1.2

Should maths be taught with many application examples? Or abstractly?

5

1.3

Notation: Some seeming trivialities . . . . . . . . . . . . . . . . . . . .

6

2 Functions & Derivatives
2.1

Basics on uni-variate functions . . . . . . . . . . . . . . . . . . . . . .

88

2.1.1 Continuous, differentiable, & smooth functions; 2.1.2 Polynomials,
piece-wise, basis functions, splines

2.2

Partial vs. total derivative, and chain rules . . . . . . . . . . . . . . . .

92.2.1 Partial Derivative; 2.2.2 Total derivative, computation graphs, forward
and backward chain rules

2.3

Gradient, Jacobian, Hessian, Taylor Expansion . . . . . . . . . . . . . .

12

2.3.1 Gradient & Jacobian; 2.3.2 Hessian; 2.3.3 Taylor expansion

2.4

Derivatives with matrices . . . . . . . . . . . . . . . . . . . . . . . . .

15

2.4.1 Derivative Rules; 2.4.2 Example: GP regression; 2.4.3 Example: Logistic regression

2.5

Check your gradients numerically! . . . . . . . . . . . . . . . . . . . .

18

2

Maths for Intelligent Systems, Marc Toussaint2.6

Further Examples and Exercises . . . . . . . . . . . . . . . . . . . . . .

19

2.6.1 More derivatives

2.7

Exercise 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

19

2.7.1 Multivariate Calculus; 2.7.2 Finite Difference Gradient Checking; 2.7.3 Backprop
in a Neural Net

2.8

Exercise 5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

20

2.8.1 Backprop in a Neural Net; 2.8.2 Logistic Regression Gradient & Hessian

3 Linear Algebra
3.1

Vector Spaces . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

21
21

3.1.1 Why should we care for vector spaces in intelligent systems research?;
3.1.2 What is a vector?; 3.1.3 What is a vector space?

3.2

Vectors, dual vectors, coordinates, matrices, tensors . . . . . . . . . . .

23

3.2.1 A taxonomy of linear functions; 3.2.2 Bases and coordinates; 3.2.3 The
dual vector space – and its coordinates; 3.2.4 Coordinates for every linear
thing: tensors; 3.2.5 Finally: Matrices; 3.2.6 Coordinate transformations

3.3

Scalar product and orthonormal basis . . . . . . . . . . . . . . . . . . .

28

3.3.1 Properties of orthonormal bases

3.4

The Structure of Transforms & Singular Value Decomposition . . . . .

31

3.4.1 The Singular Value Decomposition Theorem

3.5

Point of departure from the coordinate-free notation . . . . . . . . . .

33

3.6

Filling SVD with life . . . . . . . . . . . . . . . . . . . . . . . . . . . .

34

3.6.1 Understand vv⊤ as a projection; 3.6.2 SVD for symmetric matrices;
3.6.3 SVD for general matrices

3.7

Eigendecomposition . . . . . . . . . . . . . . . . . . . . . . . . . . . .

37

3.7.1 Power Method; 3.7.2 Power Method including the smallest eigenvalue; 3.7.3 Why should I care about Eigenvalues and Eigenvectors?

3.8

Beyond this script: Numerics to compute these things . . . . . . . . . .

38

3.9

Derivatives as 1-forms, steepest descent, and the covariant gradient . .

39

3.9.1 The coordinate-free view: A derivative takes a change-of-input vector
as input, and returns a change of output; 3.9.2 Contra- and co-variance;
3.9.3 Steepest descent and the covariant gradient vector

Maths for Intelligent Systems, Marc Toussaint3

3.10 Examples and Exercises . . . . . . . . . . . . . . . . . . . . . . . . . .

42

3.11 Exercise 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

42

3.11.1 Basis; 3.11.2 From the Robotics Course; 3.11.3 Bases for Polynomials

3.12 Exercise 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

44

3.12.1 Projections; 3.12.2 SVD; 3.12.3 Bonus: Scalar product and Orthogonality

3.13 Exercise 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

45

3.13.1 Eigenvectors; 3.13.2 Covariance and PCA; 3.13.3 Bonus: RKHS

4 Optimization
4.1

Downhill algorithms for unconstrained optimization . . . . . . . . . . .

47
47

4.1.1 Why you shouldn’t trust the magnitude of the gradient; 4.1.2 Ensuring
monotone and sufficient decrease: Backtracking line search, Wolfe conditions,& convergence; 4.1.3 The Newton direction; 4.1.4 Least Squares & GaussNewton: a very important special case; 4.1.5 Quasi-Newton & BFGS: approximating the hessian from gradient observations; 4.1.6 Conjugate Gradient;4.1.7 Rprop*

4.2

The general optimization problem – a mathematical program . . . . . .

58

4.3

The KKT conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . .

59

4.4

Unconstrained problems to tackle a constrained problem . . . . . . . .

60

4.4.1 Augmented Lagrangian*

4.5

The Lagrangian . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

62

4.5.1 How the Lagrangian relates to the KKT conditions; 4.5.2 Solving
mathematical programs analytically, on paper.; 4.5.3 Solving the dual problem, instead of the primal.; 4.5.4 Finding the “saddle point” directly with a
primal-dual Newton method.; 4.5.5 Log Barriers and the Lagrangian

4.6

Convex Problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

67

4.6.1 Convex sets, functions, problems; 4.6.2 Linear and quadratic programs; 4.6.3 The Simplex Algorithm; 4.6.4 Sequential Quadratic Programming

4.7Blackbox & Global Optimization: It’s all about learning . . . . . . . . .

71

4.7.1 A sequential decision problem formulation; 4.7.2 Acquisition Functions for Bayesian Global Optimization*; 4.7.3 Classical model-based blackbox optimization (non-global)*; 4.7.4 Evolutionary Algorithms*

4.8

Examples and Exercises . . . . . . . . . . . . . . . . . . . . . . . . . .

75

4

Maths for Intelligent Systems, Marc Toussaint4.9

Exercise 6 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

75

4.9.1 Convergence proof; 4.9.2 Backtracking Line Search; 4.9.3 GaussNewton; 4.9.4 Robust unconstrained optimization

4.10 Exercise 7 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

78

4.10.1 Lagrangian Method of Multipliers; 4.10.2 Equality Constraint Penalties
and Augmented Lagrangian

4.11 Exercise 8 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

79

4.11.1 Lagrangian and dual function; 4.11.2 Optimize a constrained problem

4.12 Exercise 9 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

80

4.12.1 Network flow problem; 4.12.2 Minimum fuel optimal control; 4.12.3 Reformulating
an ℓ1 -norm

4.13 Exercise 10 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

81

4.13.1 Restarts of Local Optima; 4.13.2 GP-UCB Bayesian Optimization

5 Probabilities & Information
5.1

Basics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

82
82

5.1.1 Axioms, definitions, Bayes rule; 5.1.2 Standard discrete distributions;5.1.3 Conjugate distributions; 5.1.4 Distributions over continuous domain;
5.1.5 Gaussian; 5.1.6 “Particle distribution”

5.2

Between probabilities and optimization: neg-log-probabilities, exp-negenergies, exponential family, Gibbs and Boltzmann . . . . . . . . . . . .

88

5.3

Information, Entropie & Kullback-Leibler . . . . . . . . . . . . . . . . .

91

5.4

The Laplace approximation: A 2nd-order Taylor of log p . . . . . . . . .

92

5.5

Variational Inference . . . . . . . . . . . . . . . . . . . . . . . . . . . .

93

5.6

The Fisher information metric: 2nd-order Taylor of the KLD . . . . . .

93

5.7

Examples and Exercises . . . . . . . . . . . . . . . . . . . . . . . . . .

94

5.8

Exercise 10 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

94

5.8.1 Maximum Entropy and Maximum Likelihood

5.9

Exercise 11 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

95

5.9.1 Maximum likelihood and KL-divergence; 5.9.2 Laplace Approximation; 5.9.3 Learning = Compression; 5.9.4 A gzip experiment; 5.9.5 Maximum
Entropy and ML

A Gaussian identities

97Maths for Intelligent Systems, Marc Toussaint5

B Further

102

Index

103

1

Speaking Maths

1.1

Describing systems

Systems can be described in many ways. Biologists describe their systems often using
text, and lots and lots of data. Architects describe buildings using drawings. Physicists
describe nature using differential equations, or optimality principles, or differential geometry and group theory. The whole point of science is to find descriptions of systems—in
the natural science descriptions that allow prediction, in the engineering sciences descriptions that enable the design of good systems, problem-solving systems.
And how should we describe intelligent systems? Robots, perception systems, machinelearning systems? I think there are two main categories: the imperative way in terms ofliteral algorithms (code), or the declarative way in terms of formulating the problem. I
prefer the latter.
The point of this lecture is to teach you to speak maths, to use maths to describe
systems or problems. I feel that most maths courses rather teach to consume maths, or
solve mathematical problems, or prove things. Clearly, this is also important. But for the
purpose of intelligent systems research, it is essential to be skilled in expressing problems
mathematically, before even thinking about solving them and deriving algorithms.
If you happen to attend a Machine Learning or Robotics course you’ll see that every
problem is addressed the same way: You have an “intuitively formulated” problem; thefirst step is to find a mathematical formulation; the second step to solve it. The second
step is often technical. The first step is really the interesting and creative part. This iswhere you have to nail down the problem, i.e., nail down what it means to be successful
or performant – and thereby describe “intelligence”, or at least a tiny aspect of it.
The “Maths for Intelligent Systems” course will recap essentials of multi-variate functions, linear algebra, optimization, and probabilities. These fields are essential to formulate problems in intelligent systems research and hopefully will equip you with the basics
of speaking maths.

1.2

Should maths be taught with many application examples? Orabstractly?

Maybe this is the wrong question and implies a view on maths I don’t agree with. I think(but this is arguable) maths is nothing but abstractions of real-world things. At least I
aim to teach maths as abstractions of real-world things. It is misleading to think that

6

Maths for Intelligent Systems, Marc Toussaintthere is “pure maths” and then “applications”. Instead mathematical concepts, such asa vector, are abstractions of real-world things, such as faces, scenes, images, documents;
and theorems, methods and algorithms that apply on vectors of course also apply to all
the real-world things—subject to the limitations of this abstraction. So, the goal is not
to teach you a lookup table of which method can be used in which application, butrather to teach which concepts maths offers to abstract real-world things—so that youfind such abstractions yourself once you’ll have to solve a real-world problem.
But yes, I believe that maths – in our context – should ideally be taught with many
exercises relating to AI problems. Perhaps the ideal would be: Teach Maths using AI exercises (where AI problems are formulated and treatedanalytically).
 Teach AI using coding exercises.
 Teach coding using maths-implementation exercises.

But I haven’t yet adopted this myself in my teaching.

1.3

Notation: Some seeming trivialities

Equations and mathematical expressions have a syntax. This is hardly ever made explicit1
and might seem trivial. But it is surprising how buggy mathematical statements can bein scientific papers (and oral exams). I don’t want to write much text about this, justsome bullet points:
 Always declare mathematical objects.
 Be aware of variable and index scoping. For instance, if you have an equation,and one side includes a variable i, the other side doesn’t, this often is a notationalbug. (Unless this equation actually makes a statement about independence on i.)
 Type checking. Within an equation, be sure to know exactly of what type eachterm is: vector? matrix? scalar? tensor? Is the type and dimension of both sidesof the equation consistent?
 Decorations are ok, but really not necessary. It is much more important to declare−v , |v⟩
all things. E.g., there are all kinds of decorations used for vectors, v, v, →and matrices. But these are not necessary. Properly declaring all symbols is much
more important.
 When declaring sets of indexed elements, I use the notation {xi }ni=1 . Similarly fortuples: (xi )ni=1 , (x1 , .., xn ), x1:n .
1 Except perhaps by Gödel’s incompleteness theorems and areas like automated theorem proving.

Maths for Intelligent Systems, Marc Toussaint7

 When defining sets, we write something like {f (x) : x ∈ R}, or {n ∈ N :∃ {vi }ni=1 linearly independent, vi ∈ V } I usually use brackets [a = b] ∈ {0, 1} for the boolean indicator function of some
expression. An alternative notation is I(a = b), or the Kronecker symbol δab .
 A tuple (a, b) ∈ A × B is an element of the product space.
 direct sum A ⊕ B is same as A × B in finite-dimensional spaces
 If f : X → Y , then minx f (x) ∈ Y is minimal function value (output); whereas
argminx f (x) ∈ X is the input (“argument”) that minimizes the function. E.g.,
minx f (x) = f (argminx f (x)).
 One should distinguish between the infimum inf x f (x) and supremum supx f (x)from the min/max: the inf/sup refer to limits, while the min/max to values actuallyacquired by the function. I must admit I am sloppy in this regard and usually only
write min/max.
 Never use multiple letters for one thing. E.g. length = 3 means l times e times ntimes g times t times h equals 3.
 There is a difference between → and 7→:

f : R → R, x 7→ cos(x)(1)

 The dot is used to help defining functions with only some arguments fixed:

f : A×B →C ,f (a, ·) : B → C, b 7→ f (a, b)(2)Another example is the typical declaration of an inner product: ⟨·, ·⟩ : V × V → R.
P The p-norm or Lp -norm is defined as ||x||p = [ i xpi ]1/p . By default p = 2, that2
is ||x|| = ||x||2 = |x|, which
P is the L -norm or Euclidean length of a vector. Also
1
the L -norm ||x||1 = i |xi | (aka Manhattan distance) plays an important role.
 I use diag(a1 , .., an ) =

a1









0


..








.

for a diagonal matrix. And overload this0anso that diag(A) = (A11 , .., Ann ) is the diagonal vector of the matrix A ∈ Rn×n .
 A typical convention is

0n = (0, .., 0) ∈ Rn ,1n = (1, .., 1) ∈ Rn ,In = diag(1n )(3)nAlso, ei = (0, .., 0, 1, 0, .., 0) ∈ R often denotes the ith column of the identity
matrix, which of course are the coordinates of a basis vector ei ∈ V in a basis(ei )ni=1 .
 The element-wise product of two matrices A and B is also called Hadamard
product and notated A ◦ B (which has nothing to do with the concatenation of
two operations). If there is need, perhaps also use this to notate the element-wiseproduct of two vectors.

8

Maths for Intelligent Systems, Marc Toussaint2

Functions & Derivatives

2.1

Basics on uni-variate functions

We super quickly recap the basics of functions R → R, which the reader might already
know.
2.1.1

Continuous, differentiable, & smooth functions

 A function f : R → R is continuous at x if the limit limh→0 f (x + h) = f (x)exists and equals f (x) (from both sides).
(x)
 A function f : R → R is differentiable at x if f ′ (x) = limh→0 f (x+h)−fexists.
hNote, differentiable ⇒ continuous.

 A function is continuous/differentiable if it is continuous/differentiable at anyx ∈ R.
 A function f is an element of Ck if it is k-fold continuously differentiable, i.e., if
its k-th derivative f (k) is continuous. For example, C0 is the space of continuousfunctions, and C2 the space of twice continuously differentiable functions.
 A function f ∈ C∞ is called smooth (infinitely often differentiable).

2.1.2

Polynomials, piece-wise, basis functions, splines

Let’s recap some basic functions that often appear in AI research.
Pp
 A polynomial of degree p is of the form f (x) = i=0 ai xi , which is a weighed
sum of monomials 1, x, x2 , ... Note that for multi-variate functions, the numberof monomials grows combinatorially with the degree and dimension. E.g., the222
monomials
 of degree 2 in 3D space are x1 , x1 x2 , x1 x3 , x2 , x2 x3 , x3 . In general,
we have dp monomials of degree p in d-dimensional space.
 Polynomials are smooth.
 I assume it is clear what piece-wise means (we have different polynomials in disjoint
intervals covering the input domain R).
P A set of basis functions {e1 , ..., en } defines the space f = { i ai ei : ai ∈ R}
of functions that are linear combinations of the basis functions. More details on
bases are discussed below for vector spaces in general. Here we note typical basis
functions:– Monomials, which form the basis of polynomials.

Maths for Intelligent Systems, Marc Toussaint9– Trigonometric functions sin(2πax), cos(2πax), with integers a, which form
the Fourier basis of functions.
– Radial basis functions φ(|x−ci |), where we have n centers ci ∈ R, i = 1, .., n,and φ is typically some bell-shaped or local support function around zero. A
very common one is the squared exponential φ(d) = exp{−ld2 } (which you
might call a non-normalized Gaussian with variance 1/l).
 The above functions are all examples of parameteric functions, which meansthat they can be specified by a finite number of parameters ai . E.g., when we
have a finite set of basis functions, the functions can all be described by the finite
set of weights in the linear combination.

However, in general a function f : R → R is an “infinite-dimensional object”,
i.e., it has infinitely many degrees-of-freedom f (x), i.e., values at infinitely many
points x. In fact, sometimes it is useful to think of f as a “vector” of elements fx with continuous index x. Therefore, the space of all possible functions,and also the space of all continuous function C0 and smooth functions C∞ , isinfinite-dimensional. General functions cannot be specified by a finite number ofparameters, and they are called non-parameteric.
The core example are functions used for regression or classification in Machine
Learning, which are a linear combination of an infinite set of basis functions. E.g.,
an infinite set of radial basis functions φ(|x − c|) for all centers c ∈ R. This
infinite set of basis functions spans a function space called Hilbert space (in the
ML context, “Reproducing Kernel Hilbert Space (RKHS)”), which is an infinitedimensional vector space. Elements in that space are called non-parameteric.
 As a final note, splines are parameteric functions that are often used in robotics
and engineering in general. Splines usually are piece-wise polynomials that arecontinuously joined. Namely, a spline of degree p is in Cp−1 , i.e., p − 1-fold continuously differentiable. A spline is not fully smooth, as the p-th derivative isdiscontinuous. E.g., a cubic spline has a piece-wise constant (“bang-bang”) jerk
(3rd derivative). Cubic splines and B-splines (which are also piece-wise polynomials) are commonly used in robotics to describe motions. In computational design
and graphics, B-splines are equally common to describe surfaces (in the form of
Non-uniform rational basis spline, NURBS). This link contain a brief reference forsplines.

2.2

Partial vs. total derivative, and chain rules

2.2.1

Partial Derivative

A multi-variate function f : Rn → R can be thought of as a function of n arguments,
f (x1 , .., xn ).

10Maths for Intelligent Systems, Marc ToussaintDefinition 2.1. The partial derivative of a function of multiple arguments f (x1 , .., xn )is the standard derivative w.r.t. only one of its arguments,
f (x1 ,..,xi +h,..,xn )−f (x)∂.
∂xi f (x1 , .., xn ) = limh→0
h

2.2.2

(4)

Total derivative, computation graphs, forward and backward chain rules

Let me start with an example: We have three real-valued quantities x, g and f which
depend on each other. Specifically,f (x, g) = 3x + 2g

and g(x) = 2x .

(5)

Question: What is the “derivative of f w.r.t. x”?
The correct answer is: Which one do you mean? The partial or total?
The partial derivative defined above really thinks of f (x, g) as a function of two arguments, and does not at all care about whether there might be dependencies of thesearguments. It only looks at f (x, g) alone and takes the partial derivative (=derivative
w.r.t. one function argument):∂∂x f (x, g) = 3(6)

However, if you suddenly talk about h(x) = f (x, g(x)) as a function of the argument x
only, that’s a totally different story, and∂∂∂x h(x) = ∂x 3x + 2(2x) = 7(7)

Bottom line, the definition of the partial derivative really depends on what you explicitly
defined as the arguments of the function.
To allow for a general treatment of differentiation with dependences we need to definea very useful concept:
Definition 2.2. A function network or computation graph is a directed acyclic
graph (DAG) of n quantities xi where each quantity is a deterministic function of
a set of parents π(i) ⊂ {1, .., n}, that isxi = fi (xπ(i) )(8)where xπ(i) = (xj )j∈π(i) is the tuple of parent values. This could also be called adeterministic Bayes net.
In a function network all values can be computed deterministically if the input values(which do have no parents) are given. Concerning differentiation, we may now ask:
Assume we have a variation dx of some input value, how do all other values vary? The
chain rules give the answer. It turns our there are two chain rules in function networks:

Maths for Intelligent Systems, Marc Toussaint

11

f

f

g
g
x

x

Figure 1:
General Chain Rule. Left: Forward-Version, Right: Backward-Version.
They gray arc denotes the direct dependence ∂f∂x , which appears in the summations
via dx/dx ≡ 1, df /df ≡ 1.

Identities 2.1 (Chain rule (Forward-Version)).
X ∂f dg
df
=
dx
∂g dx

(with

g∈π(f )

dx
≡ 1, in case x ∈ π(f ))dx(9)Read this as follows: “The change of f with x is the sum of changes that come fromits direct dependence on g ∈ π(f ), each multiplied the change of g with x.”
dfThis rule defines the total derivative of dx
w.r.t. x. Note how different these two
notions of derivatives are by definition: a partial derivative only looks at a function itself
and takes a limit of differences w.r.t. one argument—no notion of further dependencies.
The total derivative asks how, in a function network, one value changes with a changeof another.

The second version of the chain rule is:Identities 2.2 (Chain rule (Backward-Version)).
df
=
dx

X
g:x∈π(g)

df ∂g
dg ∂x

(with

df
≡ 1, in case x ∈ π(f ))
df

(10)

Read this as follows: “The change of f with x is the sum of changes that arise from allchanges of g which directly depend on x.”
Figure 1 illustrates the fwd and bwd versions of the chain rule. The bwd version allows
df
you to propagate back, given gradients dg
from top to g, one step further down, fromdg
top to x. The fwd version allows you to propagate forward, given gradients dx
from g tobottom, one step further up, from f to bottom. Both versions are recursive equations.

12

Maths for Intelligent Systems, Marc ToussaintIf you would recursively plug in the definition for a given functional network, both ofdfthem would yield the same expression of dx
in terms of partial derivatives only.
Let’s compare to the chain rule as it is commonly found in other texts (written more
precisely):∂f (g)∂g(x)∂f (g(x))=∂x∂g g=g(x) ∂x(11)Note that we here very explicitly notated that ∂f∂g(g) considers f to be a function of
the argument g, which is evaluted at g = g(x). Written like this, the rule is fine. But
the above discussion and explicitly distinguishing between partial and total derivative is,when things get complicated, less prone to confusion.

2.3

Gradient, Jacobian, Hessian, Taylor Expansion

2.3.1

Gradient & JacobianLet’s take the next step and consider functions f : Rn → Rd that map from n numbers
to a d-dimensional output. In this case, we can take the partial derivative of each output
w.r.t. each input argument, leading to a matrix of partial derivatives:
Definition 2.3. Given f : Rn → Rd , we define the derivative (also called Jacobianmatrix) as∂

f1 (x)
f2 (x)∂..
∂x f (x) =.
f∂x1 d (x)


 ∂x1


 ∂

 ∂x1






 ∂

∂
∂x2 f1 (x)∂∂x2 f2 (x)

...
...

∂
∂x2 fd (x)

...


∂
∂xn f1 (x) ∂∂xn f2 (x) 




..



.


∂
f
(x)
d
∂xn

(12)

When the function only as one output dimension, f : Rn → R1 , the partial derivative
can be written as a vector. Unlike many other texts, I advocate for consistency with theJacobian matrix (and contra-variance, see below) and define this to be a row vector:
Definition 2.4. We define the derivative of f : Rn → R as the row vector ofpartial derivatives∂∂∂x f (x) = ( ∂x1 f,.., ∂x∂n f ) .

(13)

Further, we define the gradient as the corresponding column “vector”∇f (x) =hi⊤∂∂x f (x).

(14)

Maths for Intelligent Systems, Marc Toussaint13

The “purpose” of a derivative is to output a change of function value when being
multiplied to a change of input δ. That is, in first order approximation, we havef (x + δ) − f (x) =˙ ∂f (x) δ ,(15)

where =
˙ denotes “in first order approximation”. This equation holds, no matter if theoutput space is Rd or R, or the input space and variation is δ ∈ Rn or δ ∈ R. In thegradient notation we havef (x + δ) − f (x) =˙ ∇f (x)⊤δ .

(16)

Jumping ahead to our later discussion of linear algebra: The above two equations are
written in coordinates. But note that the equations are truly independent of the choiceof vector space basis and independent of an optional metric or scalar product in V . The
transpose should not be understood as a scalar product between two vectors, but rather
as undoing the transpose in the definition of ∇f . All this is consistent to understanding
the derivatives as coordinates of a 1-form, as we will introduce it later.
Given a certain direction d (with |d| = 1) we define the directional derivative as∇f (x)⊤d, and it holdsf (x + ϵd) − f (x).
ϵ→0ϵ∇f (x)⊤d = lim

2.3.2

(17)

HessianDefinition 2.5. We define the Hessian of a scalar function f : Rn → R as thesymmetric matrix
∂∇2 f (x) = ∂x
∇f (x) =∂2

 ∂x ∂x

1
1


∂2

 ∂x ∂x
2
1







2


f
f

..
.

∂
∂xn ∂x1 f∂2∂x1 ∂x2 f∂2∂x2 ∂x2 f∂2∂xn ∂x2 f

...
...

∂2
∂x1 ∂xn f 



∂2

f
∂x2 ∂xn 



..
.

...

∂2
∂xn ∂xn f









(18)

The Hessian can be thought of as the Jacobian of ∇f . Using the Hessian, we can express
the 2nd order approximation of f as:
1
f (x + δ) =¨ f (x) + ∂f (x) δ + δ⊤ ∇2 f (x) δ .
2(19)

For a uni-variate function f : R → R, the Hessian is just a single number, namely thesecond derivative f ′′ (x). In this section, let’s call this the “curvature” of the function(not to be confused with the Rimannian curvature of a manifold). In the uni-variatecase, we have the obvious cases:

14

Maths for Intelligent Systems, Marc Toussaint If f ′′ (x) > 0, the function is locally “curved upwared” and convex (see also theformal Definition 4.6.1).
 If f ′′ (x) < 0, the function is locally “curved downward” and concave.

n
In the multi-variate caseP f : R → R, the Hessian matrix H is symmetric and we candecompose it as H = i λi hi h⊤i with eigenvalues λi and eigenvectors hi (which we willlearn about in detail later). Importantly, all hi will be orthogonal to each other, forminga nice orthonormal basis.

This insight gives us a very strong intuition on how the Hessian H describes the local
curvature of the function f : λi gives the directional curvature, i.e., the curvature in
the direction of eigenvector hi . If λi > 0, f is curved upward along hi ; if λi < 0, f iscurved downward along hi . Therefore, the eigenvalues λi tell us whether the function is
locally curved upward, downward, or flat in each of the orthogonal directions hi .
∂f = 0 is zero, i.e., the derivative (slope) of theThis becomes particular intuitive if ∂xfunction is zero in all directions. When the curvatures λi are positive in all directions,the function is locally convex (upward parabolic) and x is a local minimum; if thecurvatures λi are all negative, the function is concave (downward parabolic) and x is alocal maximum; if some curvatures are positive and some are negative along different
directions hi , then the function curves down in some directions, and up in others, andx is a saddle point.

Again jumping ahead, in the coordinate-free notation, the second derivative would be
defined as the 2-form
d2 f x : V × V → G,(v, w) 7→ lim

df x+hw (v) − f x (v)hf (x + hw + lv) − f (x + hw) − f (x + lv) + f (x)
= lim
.
h,l→0
hl
h→0

(20)
(21)
(22)

The Hessian matrix are the coordinates of this 2-form (which would actually be a row
vector of row vectors).
2.3.3

Taylor expansion

In 1D, we have11
f (x + v) ≈ f (x) + f ′ (x)v + f ′′ (x)v 2 + · · · + f (k) (x)v k
2
k!

(23)

For f : Rn → R, we have1f (x + v) ≈ f (x) + ∇f (x)⊤v + v⊤∇2 f (x)v + · · ·2(24)Maths for Intelligent Systems, Marc Toussaint15

which is equivalent tof (x + v) ≈ f (x) +X1∂∂xj f (x)vj + 2∂2
jk ∂xj ∂xk f (x)vj vk + · · ·

P

(25)

j

2.4

Derivatives with matrices

The next section will introduce linear algebra from scratch – here we first want to learn
how to practically deal with derivatives in matrix expressions. We think of matrices andvectors simply as arrays of numbers ∈ Rn×m and Rn . As a warmup, try to solve thefollowing exercises:(i) Let X, A be arbitrary matrices, A invertible. Solve for X:
XA + A⊤ = I(26)(ii) Let X, A, B be arbitrary matrices, (C − 2A⊤) invertible. Solve for X:
X⊤C = [2A(X + B)]⊤

(27)(iii) Let x ∈ Rn , y ∈ Rd , A ∈ Rd×n . A obviously not invertible, but let A⊤A be
invertible. Solve for x:(Ax − y)⊤A = 0⊤n(28)(iv) As above, additionally B ∈ Rn×n , B positive-definite. Solve for x:(Ax − y)⊤A + x⊤B = 0⊤n(29)(v) A core problem in Machine Learning: For β ∈ Rd , y ∈ Rn , X ∈ Rn×d , computeargmin ||y − Xβ||2 + λ||β||2 .

(30)

β

(vi) A core problem in Robotics: For q, q0 ∈ Rn , ϕ : Rn → Rd , y ∗ ∈ Rd non-linear
but smooth, computeargmin ||ϕ(q) − y ∗ ||2C + ||q − q0 ||2W .

(31)

q

Use a local linearization of ϕ to solve this.
For problem (v), we want to find a minimum for a matrix expression. We find this bysetting the derivative equal to zero. Here is the solution, and below details will become
clear:∂
0 = ∂β||y − Xβ||2 + λ||β||2

(32)

16

Maths for Intelligent Systems, Marc Toussaint= 2(y − Xβ)⊤(−X) + 2λβ⊤

(33)

⊤

0 = −X (y − Xβ) + λβ
⊤

(34)

⊤

0 = −X y + (X X + λI)β(35)

β = −(X⊤X + λI)-1 X⊤y(36)Line 2 uses a standard rule for the derivative (see below) and gives a row vector equation.
Line 3 transposes this to become a column vector equation.
2.4.1

Derivative Rules

As 2nd order terms are very common in AI methods, this is a very useful identity tolearn:
Identities 2.3.
⊤ ∂
∂
⊤
⊤ ⊤∂
∂x f (x) Ag(x) = f (x) A ∂x g(x) + g(x) A ∂x f (x)(37)

Note that using the ’gradient column’ convention this reads∂∂g(x)]⊤A⊤f (x) + [ ∂xf (x)]⊤Ag(x)∇x f (x)⊤Ag(x) = [ ∂x(38)which I find impossible to remember, and mixes gradients-in-columns (∇) with gradientsin-rows (the Jacobian) notation.
Special cases and variants of this identity are:∂∂x [whatever]x = [whatever] ,⊤∂ ⊤∂x a x = a∂∂x Ax = A⊤⊤

if whatever is indep. of x⊤∂∂x (Ax − b) (Cx − d) = (Ax − b) C + (Cx − d) A
⊤
⊤ ⊤
∂ ⊤
∂x x Ax = x A + x A∂∂
1 ⊤ − 121 ⊤⊤ 12
2x⊤ = ||x||x∂x ||x|| = ∂x (x x) = 2 (x x)
2

⊤
⊤ ⊤
∂
⊤
∂x2 (Ax + a) C(Bx + b) = A CB + B C A

(39)
(40)
(41)
(42)
(43)
(44)
(45)

Further useful identities are:Identities 2.4 (Derivative Rules).
-1 ∂∂∂θ |A| = |A| tr(A
∂θ A)
-1
∂
∂
-1
-1
∂θ A = −A ( ∂θ A) A

(46)
(47)

Maths for Intelligent Systems, Marc Toussaint∂∂θ tr(A) =

X

17

∂
∂θ Aii(48)

i

We can also directly take a derivative of a scalar value w.r.t. a matrix:⊤∂ ⊤∂X a Xb = ab
⊤
⊤ ⊤

(49)

⊤
⊤
∂
∂X (a X CXb) = C Xab + CXba∂
∂X tr(X) = I

(50)
(51)

But if this leads to confusion I would recommend to never take a derivative w.r.t. amatrix. Instead, perhaps take the derivative w.r.t. a matrix element: ∂X∂ ij a⊤Xb = ai bj .
For completeness, here are the most important matrix identities (the appendix lists
more):
Identities 2.5 (Matrix Identities).
(A-1 + B -1 )-1 = A (A + B)-1 B = B (A + B)-1 A
-1

-1 -1

(A − B ) = A (B − A)

-1

B

(53)

(A + U BV )-1 = A-1 − A-1 U (B -1 + V A-1 U )-1 V A-1
-1

-1 -1

-1

(A + B ) = A − A(B + A) A
⊤

-1 ⊤

⊤

-1

-1 ⊤

-1

(54)
(55)

-1 ⊤ -1(A + J BJ) J B = A J (B + JA J )
⊤

(52)

-1 ⊤

(A + J BJ) A = I − (A + J BJ) J BJ

(56)
(57)

(54)=Woodbury; (56,57) holds for pos def A and B. See also the matrix cookbook.

2.4.2

Example: GP regression

An example from GP regression: The log-likelihood gradient w.r.t. a kernel hyperparameter:1n1log P (y|X, b) = − y⊤K -1 y − log |K| − log 2π2222where Kij = e−b(xi −xj ) + σ 2 δij⊤-1 ∂∂ ⊤ -1-1⊤-1-1∂b y K y = y (−K ( ∂b K)K ) = y (−K AK ) ,2
with Aij = −(xi − xj )2 e−b(xi −xj )

1 ∂
∂
1
-1 ∂
-1
∂b log |K| = |K| ∂b |K| = |K| |K| tr(K ∂b K) = tr(K A)

(58)
(59)

(60)
(61)

18

Maths for Intelligent Systems, Marc Toussaint2.4.3

Example: Logistic regression

An example from logistic regression: We have the loss gradient and want the Hessian:∇β L = X⊤(p − x) + 2λIβ
wherepi = σ(x⊤i β) ,(62)z

σ(z) =e,
1 + ezσ ′ (z) = σ(z) (1 − σ(z))

(63)

∂
∂
∇β2 L = ∂β∇β L = X⊤ ∂β
p + 2λ(64)⊤∂∂β pi = pi (1 − pi ) xin∂∂β p = diag([pi (1 − pi )]i=1 ) X = diag(p ◦ (1 − p)) X(65)∇β2 L = X⊤diag(p ◦ (1 − p)) X + 2λI

(67)

(66)

(Where ◦ is the element-wise product.)

2.5

Check your gradients numerically!

This is your typical work procedure when implementing a Machine Learning or AI’ish or
Optimization kind of methods:
 You first mathematically (on paper/LaTeX) formalize the problem domain, including the objective function.
 You derive analytically (on paper) the gradients/Hessian of your objective function.
 You implement the objective function and these analytic gradient equations in
Matlab/Python/C++, using linear algebra packages.
 You test the implemented gradient equations by comparing them to a finite
difference estimate of the gradients!
 Only if that works, you put everything together, interfacing the objective & gradient equations with some optimization algorithm

Algorithm 1 Finite Difference Jacobian Check
Input: x ∈ Rn , function f : Rn → Rn , function df : Rn → Rd×nˆ ∈ Rd×m , and ϵ = 10−6
1: initialize J
2: for i = 1 : n do
3:Jˆ·i = [f (x + ϵei ) − f (x − ϵei )]/2ϵ// assigns the ith column of Jˆ
4: end forˆ − df (x)||∞ < 10−4 return true; else false
5: if ||J
Here ei is the ith standard basis vector in Rn .

Maths for Intelligent Systems, Marc Toussaint

19

2.6

Further Examples and Exercises

2.6.1

More derivativesa) In 3D, note that a × b = skew(a)b = −skew(b)a, where skew(v) is the skew matrixof v. What is the gradient of (a × b)2 w.r.t. a and b?

2.7

Exercise 4

2.7.1

Multivariate Calculus

Given tensors y ∈ Ra×...×z and x ∈ Rα×...×ω where y is a function of x, the Jacobian
tensor J = ∂x y is in Ra×...×z×α×...×ω and has coefficients
Ji,j,k,...,l,m,n... =∂∂xl,m,n,...

yi,j,k,...

(All “output” indices come before all “input” indices.)
Compute the following Jacobian tensors∂
x, where x is a vector(i) ∂x∂ ⊤(ii) ∂x
x Ax, where A is a matrix∂ ⊤y Ax, where x and y are vectors(iii) ∂A

(note, we take the derivative w.r.t. A)∂(iv) ∂A
Ax∂(v) ∂xf (x)⊤Ag(x), where f and g are vector-valued functions

2.7.2

Finite Difference Gradient Checking

The following exercises will require you to code basic functions and derivatives. You cancode in your prefereed language (Matlab, NumPy, Julia, whatever).
(i) Implement the following pseudo code for empirical gradient checking in the programming language of your choice:Input: x ∈ Rn , function f : Rn → Rd , function df : Rn → Rd×n
1: initialize Jˆ ∈ Rd×n , and ϵ = 10−6
2: for i = 1 : n do
3:Jˆ·i = [f (x + ϵei ) − f (x − ϵei )]/2ϵ// assigns the ith column of Jˆ
4: end for
5: if ||Jˆ − df (x)||∞ < 10−4 return true; else false

20

Maths for Intelligent Systems, Marc ToussaintHere ei is the ith standard basis vector in Rn .

(ii) Test this for f : x 7→ Ax, df : x 7→ A, where you sample x ∼randn(n,1) (N(0, 1) in
Matlab) and A ∼randn(m,n)
 f : x 7→ x⊤x, df : x 7→ 2x⊤, where x ∼randn(n,1)

2.7.3

Backprop in a Neural Net

Consider the function
f : Rh0 → Rh3 ,f (x0 ) = W2 σ(W1 σ(W0 x0 ))where Wl ∈ Rhl+1 ×hl and σ : R → R is a differentiable activation function which isapplied element-wise. We also describe the function as the computation graph:
x0 7→ z1 = W0 x0 7→ x1 = σ(z1 ) 7→ z2 = W1 x1 7→ x2 = σ(z2 ) 7→ f = W2 x2df. (Ideally also for deeper networks.)
Derive pseudo code to efficiently compute dx
0

2.8

Exercise 5

2.8.1

Backprop in a Neural Net

We consider again the function
f : Rh0 → Rh3 ,f (x0 ) = W2 σ(W1 σ(W0 x0 )) ,where Wl ∈ Rhl+1 ×hl and σ(z) = 1/(e−z + 1) is the sigmoid function which is appliedelement-wise. We established last time thatdf∂f ∂x2 ∂z2 ∂x1 ∂z1=
dx0∂x2 ∂z2 ∂x1 ∂z1 ∂x0
with:

∂xl
= diag(xl ◦ (1 − xl )) ,
∂zl

∂zl+1
= Wl ,
∂xl

∂f
= W2
∂x2

Note: In the following we still let f be a h3 -dimensional vector. For those that are
confused with the resulting tensors, simplify to f being a single scalar output.
(i) Derive also the necessary equations to get the derivative w.r.t. the weight matrices
Wl , that is the Jacobian tensor
dfdWlMaths for Intelligent Systems, Marc Toussaint

21

df
df
(ii) Write code to implement f (x) and dx
and dW.
0l

To test this, choose layer sizes (h0 , h1 , h2 , h3 ) = (2, 10, 10, 2), i.e., 2 input and 2
output dimensions, and hidden layers of dimension 10.
For testing, choose random inputs sampled from x ∼randn(2,1)
And choose random weight matrices Wl ∼ √ 1hl+1

rand(h[l+1],h[l]).

Check the implemented Jacobian by comparing to the finite difference approximation.
Debugging Tip: If your first try does not work right away, the typical approach to
debug is to “comment out” parts of your function f and df . For instance, startwith testing f (x) = W0 x0 ; then test f (x) = σ(W0 x0 ); then f (x) = W1 σ(W0 x0 );then I’m sure all bugs are found.
(iii) Bonus: Try to train the network to become the identity mapping. In the simplest case, use “stochastic gradient descent”, meaning that you sample an input,
2compute the gradients wl = d(f (x)−x), and make tiny updates Wl ← Wl − αwl .
dWl
2.8.2

Logistic Regression Gradient & Hessian

Consider the function
L : Rd → R : L(β) = −

n h
X

i
yi log σ(x⊤i β) + (1 − yi ) log[1 − σ(x⊤i β)] + λβ⊤β ,i=1where xi ∈ Rd is the ith row of a matrix X ∈ Rn×d , and y ∈ {0, 1}n is a vector of 0s and1s only. Here, σ(z) = 1/(e−z + 1) is the sigmoid function, with σ ′ (z) = σ(z)(1 − σ(z)).
∂
Derive the gradient ∂β
L(β), as well as the Hessian∇2 L(β) =∂2
L(β) .
∂β 23Linear Algebra

3.1

Vector Spaces

3.1.1

Why should we care for vector spaces in intelligent systems research?

We want to describe intelligent systems. For this we describe systems, or aspects ofsystems, as elements of a space:

22

Maths for Intelligent Systems, Marc Toussaint– The input space X, output space Y in ML– The space of functions (or classifiers) f : X → Y in ML– The space of world states S and actions A in Reinforcement Learning– The space of policies π : S → A in RL– The space of feedback controllers π : x 7→ u in robot control– The configuration space Q of a robot– The space of paths x : [0, 1] → Q in robotics– The space of image segmentations s : I → {0, 1} in computer vision
Actually, some of these spaces are not vector spaces at all. E.g. the configuration spaceof a robot might have ‘holes’, be a manifold with complex topology, or not even that(switch dimensionality at some places). But to do computations in these spaces one
always either introduces (local) parameterizations that make them a vector space,2 or
one focusses on local tangent spaces (local linearizations) of these spaces, which arevector spaces.
Perhaps the most important computation we want to do in these spaces is taking
derivatives—to set them equal to zero, or do gradient descent, or Newton steps for
optimization. But taking derivatives essentially requires the input space to (locally) be
a vector space.3 So, we also need vector spaces because we need derivatives, and Linear
Algebra to deal with the resulting equations.
3.1.2

What is a vector?

A vector is nothing but an element of a vector space.
It is in general not a column, array, or tuple of numbers. (But tuples of numbers are a
special case of a vector space.)
3.1.3

What is a vector space?

Definition 3.1 (vector space). A vector spacea V is a space (=set) on which two
operations, addition and multiplication, are defined as follows
 addition + : V × V → V is an abelian group, i.e.,– a, b ∈ V ⇒ a + b ∈ V(closed under +)– a + (b + c) = (a + b) + c (association)– a+b=b+a(commutation)

2 E.g. by definition an n-dimensional manifold X is locally isomorphic to Rn .

3 Also when the space is actually a manifold; the differential is defined as a 1-form on the local
tangent.

Maths for Intelligent Systems, Marc Toussaint– ∃ unique 0 ∈ V s.t. ∀v ∈ V : 0 + v = v

23(identity)– ∀v ∈ V : ∃ unique − v s.t. v + (−v) = 0 (inverse) multiplication · : R × V → V fulfils, for α, β ∈ R,

– α(βv) = (αβ)v– 1v = v(association)(identity)– α(v + w) = αv + αw

(distribution)a We only consider vector spaces over R.

Roughly, this definition says that a vector space is “closed under linear operations”,
meaning that we can add and scale vectors and they remain vectors.

3.2

Vectors, dual vectors, coordinates, matrices, tensors

In this section we explain what might be obvious: that once we have a basis, we canwrite vectors as (column) coordinate vectors, 1-forms as (row) coordinate vectors, andlinear transformations as matrices. Only the last subsection becomes more practical,
refers to concrete exercises, and explains how in practise not to get confused about basis
transforms and coordinate representations in different bases. So a practically oriented
reader might want to skip to the last subsection.
3.2.1

A taxonomy of linear functions

For simplicity we consider only functions involving a single vector space V . But all thatis said transfers to the case when multiple vector spaces V, W, ... were involved.
Definition 3.2. f : V → X linear ⇔ f (αv + βw) = αf (v) + βf (w), where X
is any other vector space (e.g. X = R, or X = V × V ).

Definition 3.3. f : V × V × · · · × V → X multi-linear
input.

⇔

f is linear in eachMany names are used for special linear functions—let’s make some explicit:– f : V → R, called linear functional4 , or 1-form, or dual vector.
– f : V → V , called linear function, or linear transform, or vector-valued 1-form– f : V × V → R, called bilinear functional, or 2-form– f : V × V × V → R, called 3-form (or unspecifically ’multi-linear functional’)
4 The word ’functional’ instead of ’function’ is especially used when V

is a space of functions.

24

Maths for Intelligent Systems, Marc Toussaint– f : V × V → V , called vector-valued 2-form (or unspecifically ’multi-linear map’)– f : V × V × V → V × V , called bivector-valued 3-form– f : V k → V m , called m-vector-valued k-form
This gives us a simple taxonomy of linear functions based on how many vectors a function
eats, and how many it outputs. To give examples, consider some space X of systems(examples above), which might itself not be a vector space. But locally, around a specificx ∈ X, its tangent V is a vector space. Then
– f : X → R could be a cost function over the system space.
– The differential df |x : V → R is a 1-form, telling us how f changes when ‘making
a tangent step’ v ∈ V .
– The 2nd derivative d2 f |x : V × V → R is a 2-form, telling us how df |x (v) changes
when ‘making a tangent step’ w ∈ V .
– The inner product ⟨·, ·⟩ : V × V → R is a 2-form.
Another example:– f : Ri → Ro is a neural network that maps i input signals to o output signals.
– Its derivative df |x : Ri → Ro is a vector-valued 1-form, telling us how each output
changes with a step v ∈ Ri in the input.
– Its 2nd derivative d2 f |x : Ri × Ri → Ro is a vector-valued 2-form.
This is simply to show that vector-valued functions, 1-forms, and 2-forms are common.
Instead of being a neural network, f could also be a mapping from one parameterizationof a system to another, or the mapping from the joint angles of a robot to its hand
position.
3.2.2

Bases and coordinates

We need to define some notions. I’m not commenting on these definitions—train yourself
in reading maths...
Definition 3.4. span({vi }ki=1 ) = {Pi αi vi : αi ∈ R}

Definition 3.5. {vi }ni=1 linearly independent ⇔

hP

i
α
v
=
0
⇒
∀
α
=
0
i i
i i i

Definition 3.6. dim(V ) = maxn {n ∈ N : ∃ {vi }ni=1 lin.indep., vi ∈ V }Maths for Intelligent Systems, Marc ToussaintDefinition 3.7. B = (ei )ni=1 is a basis of V25
⇔ span(B) = V and B lin.indep.

n
Definition 3.8. The tuple
P (v1 , v2 , .., vn ) ∈ R is called coordinates of v ∈ V innthe basis (ei )i=1 iff v = i vi ei

Note that Rn is also a vector space, and therefore coordinates v1:n ∈ Rn are also vectors,
but in Rn , not V . So coordinates are vectors, but vectors in general not coordinates.
P
Given a basis (ei )ni=1 , we can describe every vector v as a linear combination v = i vi ei
of basic elements—the basis vectors ei . This general idea, that “linear things” can be
described as linear combinations of “basic elements” carries over also to functions. In
fact, to all the types of functions we described above: 1-forms, 2-forms, bi-vector-valued
k-forms, whatever. And if we describe all these als linear combinations of basic elements
we automatically also introduce coordinates for these things. To get there, we first haveto introduce a second type of “basic elements”: 1-forms.
3.2.3

The dual vector space – and its coordinates

Definition 3.9. Given V , its dual space is V ∗ = {f : V → R linear} (the space of
1-forms). Every v ∗ ∈ V ∗ is called 1-form or dual vector (sometimes also covector ).
First, it is easy to see that V ∗ is also a vector space: We can add two linear functionals,
f = f1 + f2 , and scale them, and it remains a linear functional.
Second, given a basis (ei )ni=1 of V , we define a corresponding dual basis (bei )ni=1 of V ∗simply by∀i,j :ebi (ej ) = δij(68)

where δij = [i = j] is the Kronecker delta. Note that∀v ∈ V :ebi (v) = vi

(69)

That is, ebi is the 1-form that simply maps a vector to its ith coordinate. It can be shownthat (bei )ni=1 is in fact a basis of V ∗ . (Omitted.) That tells us a lot!
dim(V ∗ ) = dim(V ). That is, the space of 1-forms has the same dimension as V . At
this place, geometric intuition should kick in: indeed, every linear function over V could
be envisioned as a “plane” over V . Such a plane can be illustrated by its iso-lines andthese can be uniquely determined by their orientation and distance (same dimensionality
as V itself). Also, (assuming we’d know already what a transpose or scalar product is)
every 1-form must be of the form f (v) = c⊤v for some c ∈ V —so every f is uniquely
described by a c ∈ V . Showing that the vector space V and its dual V ∗ are really twins.
The dual basis (bei )ni=1 introduces coordinates in the dual space: Every 1-form f can be
described as a linear combination of basis 1-forms,
Xf=fi ebi(70)
i26

Maths for Intelligent Systems, Marc Toussaintwhere the tuple (f1 , f2 , .., fn ) are the coordinates of f . Andspan({bei }ni=1 ) = V ∗ .
3.2.4

(71)

Coordinates for every linear thing: tensors

We now have the basic elements: the basis vectors (ei )ni=1 of V , and basis 1-forms(bei )ni=1 of V ∗ . From these, we can describe, for instance, any bivector-valued 3-form asa linear combination as follows:f : V ×V ×V →V ×V
X
f=fijklm ei ⊗ ej ⊗ ebk ⊗ ebl ⊗ ebm

(72)
(73)

ijklm

The ⊗ is called outer product (or tensor product), and v ⊗ w ∈ V × W if V and W
are finite vector spaces. For our purposes, we may think of v ⊗ w = (v, w) simply as
the tuple of both. Therefore ei ⊗ ej ⊗ ebk ⊗ ebl ⊗ ebm is a 5-tuple and we have in total
n5 such basis objects—and fijklm denotes the corresponding n5 coordinates. The firsttwo indices are contra-variant, the last three covariant—these notions are explained indetail later.
3.2.5

Finally: Matrices

As a special case of the above, every f : V → U can be described as a linear combinationf=X

fij ei ⊗ ebj ,

(74)

ij

where (bej )nj=1 is a basis of V ∗ and (ei ) a basis of U .
Let’s see how this fits with some easier view, without all this fuss about 1-forms. We
already understood that the operator ebj (v) = vj simply picks the jth coordinate of avector. ThereforehXiXf (v) =
fij ei ⊗ ebj (v) =fij ei vj .
(75)
ij

ij

In case it helps, we can ‘derive’ this more slowly as
i
X
X
X hX
f (v) = f (vk e k ) =vk f (ek ) =vk
fij ei ⊗ ebj ek
k

=

X
ijk

k

fij vk ei ebj (ek ) =

k

X
ijk

(76)

ij

fij vk ei δjk =

XhX
i

j

i
fij vj ei .

(77)

Maths for Intelligent Systems, Marc Toussaint27

As a result, this tells us that the vector u = f (v) ∈ V has the coordinates
P ui =
And the vector f (ej ) ∈ V has the coordinates fij , that is, f (ej ) = i fij ei .

P

j fij vj .

So there are n2 coordinates fij for a linear function f . The first index is contra-variant,the second covariant (explained later). As it so happens, the whole world has agreed
on a convention on how to write such coordinate numbers on sheets of 2-dimensional
paper: as a matrix!













f11
f21
..
.
fn1

f12
f22

···
···

fn2

···

f1n 

f2n 



.. 

. 


fnn


(78)

The first (contra-variant) index spans columns; the second (covariant) spans rows. We
call this and the respective definition of a matrix multiplication as the matrix convention.
Note that the identity map I : V → V can be written as
XI=
ei ⊗ ebi , Iij = δij .

(79)

i

Equally, the (contra-variant) coordinates of a vector are written as columns













v1 

v2 



.. 

. 


vn


(80)

and the (covariant) coordinates of a 1-form h : V → R as a row
···

h2

(h1

(81)

hn )

u = f (v) is itself a vector, and its coordinates written as a column are













u1   f11


f
u2 

 21



 =  .
.. 

 .
 .
. 




un
fn1




f12
f22

···
···

fn2

···

f1n   v1 
 



f2n 
  v2    . ..   . 
 . 
. 
 

 

fnn vn
 



(82)

where this matrix multiplication is defined by ui =columns
vectoroutput spaceco-variant
contra-variant coordinatesPj fij vj , consistent to the above.

rows
1-formco-vector
derivative
input spacecontra-variant
co-variant coordinates

28

Maths for Intelligent Systems, Marc Toussaint3.2.6

Coordinate transformations

The above was rather abstract. The exercises demonstrate representing vectors andtransformations with coordinates and matrices in different input and output bases. We
just summarize here: We have two bases A = (a1 , .., an ) and B = (b1 , .., bn ), and the transformation
T that maps each ai to bi , i.e., B = T A.
 Given a vector x we denote its coordinates in A by [x]A or briefly as xA . And
we denotes its coordinates in B as [x]B or xB . E.g., xA
i is the ith coordinate inbasis A.
 [bi ]A are the coordinates of the new basis vectors in the old basis. The coordinatetransformation matrix B is given with elements Bij = [bj ]A
i . Note that[x]A = B[x]B ,(83)

i.e., while the basis transform T carries old basis ai vectors to new basis vectors bi ,the matrix B carries coordinates [x]B in the new basis to coordinates [x]A in theold basis! This is the origin of understanding that coordinates are contra-variant.
 Given a linear transform f in the vector space, we can represent it as a matrix
in four ways, using basis A or B in the input and output spaces, respectively. If[f ]AA = F is the matrix in old coordinates (using A for input and output), then[f ]BB = B -1 F B is its matrix in new coordinates, [f ]AB = F B is its matrix using
B for the input and A for the output space, and [f ]BA = B -1 F is the matrix using
A for input and B for output space.
 T itself is also a linear transform. [T ]AA = B is its matrix in old coordinates.
And the same [T ]BB = B is also its matrix in new coordinates! [T ]BA = I is its
matrix when using A for input and B for output space. And [T ]AB = B 2 is its
matrix using B for input and A for output space.

3.3

Scalar product and orthonormal basis

Please note that so far we have not in any way referred to a scalar product or a transpose.
All the concepts above, dual vector space, bases, coordinates, matrix-vector multiplication, are fully independent of the notion of a scalar product or transpose. Columns
and rows naturally appear as coordinates of vectors and 1-forms. But now we need to
introduce scalar products.
Definition 3.10. A scalar product (also called inner product) of V is a symmetricpositive definite 2-form⟨·, ·⟩ : V × V → R .

(84)

Maths for Intelligent Systems, Marc Toussaint29

with ⟨v, w⟩ = ⟨w, v⟩ and ⟨v, v⟩ > 0 for all v ̸= 0 ∈ V .
Definition 3.11. Given a scalar product, we define for every v ∈ V its dual v ∗ ∈ V ∗asX
Xv ∗ = ⟨v, ·⟩ =vi ⟨ei , ·⟩ =vi e∗i .
(85)
iiNote that ebi and e∗i are in general different 1-forms! The canonical dual basis (bei )ni=1
is independent of an introduction of a scalar product, they were the basis to introduce
coordinates for linear functions, including matrices. And while such coordinates dodepend on a choice of basis (ei )ni=1 , they do not depend on a choice of scalar product.
The 1-forms (e∗i )ni=1 also form a basis for V ∗ , but a different one to the canonicalbasis, and one that depends on the notion of a scalar product. You can see this: thecoordinates vi of v ∗ in the basis (e∗i )ni=1 are identical to the coordinates vi of v in thebasis (ei )ni=1 , but different to the coordinates (v ∗ )i of v ∗ in the basis (bei )ni=1 .
Definition 3.12. Given a scalar product, a set of vectors {vi }ni=1 is called orthonormal iff ⟨vi , vj ⟩ = δij .
Definition 3.13. Given a scalar product and basis (ei )ni=1 , we define the metrictensor gij = ⟨ei , ej ⟩, which are the coordinates of the 2-form ⟨·, ·⟩, that isX⟨·, ·⟩ =gij ebi ⊗ ebj .
(86)
ijThis also implies thatX
X
⟨v, w⟩ =vi wj ⟨ei , ej ⟩ =vi wj gij = v⊤Gw .
ij

(87)

ij

Although related, do not confuse gij with the usual definition of a metric d(·, ·) in a
metric space.
3.3.1

Properties of orthonormal bases

If we have an orthonormal basis (ei )ni=1 , many thing simplify a lot. Throughout thissubsection, we assume {ei } orthonormal.
 The metric tensor gij = ⟨ei , ej ⟩ = δij is the identity matrix.5 Such a metric isalso called Euclidean. The norm ||ϵi || = 1. The canonical dual basis (bei )ni=1 and5 Being picky, a metric is not a matrix but a twice covariant tensor (a row of rows). That’s why it iscorrectly called metric tensor.

30

Maths for Intelligent Systems, Marc Toussaintthe one defined via the scalar product (e∗i )ni=1 become identical, ebi = e∗i = ⟨ei , ·⟩.
Consequently, v and v ∗ have the same coordinates vi = (v ∗ )i w.r.t. (ei )ni=1 and(bei )ni=1 , respectively.
 The coordinates of vectors can now easily been computed:
XXXv=vi ei ⇒ ⟨ei , v⟩ = ⟨ei ,vj e j ⟩ =⟨ei , ej ⟩vj = vi
i

j

(88)

j

 The coordinates of a linear transform can equally easily been computed: Givena linear transform f : V → U an arbitrary (e.g. non-orthonormal) input basis(vi )ni=1 of V , but an orthonormal basis (ui )ni=1 , thenXf=fij ui ⊗ vbj ⇒(89)ij⟨ui , f vj ⟩ = ⟨uj ,X

fkl uk ⊗ vbl (vj )⟩ =

kl

=

X

X

fkl ⟨uj , uk ⟩ vbl (vj )klfkl δjk δlj = fij

(90)

kl

 The projection onto a basis vector is given by ei ⟨ei , ·⟩.
 The projection onto the span of several basis vectors (e1 , .., ek ) is given byPdim(V ) The identity mapping I : V → V is given by I = i=1ei ⟨ei , ·⟩.
 The scalar product with an orthonormal basis isXX⟨v, w⟩ =vi wj δij =
vi wi
ij

Pk

i=1 ei ⟨ei , ·⟩.

(91)

i

which, using matrix convention, can also be written as
w1 

X
w2 

⊤
 = v w =⟨v, w⟩ = (v1 v2 .. vn ) . 
vi w i ,


.. 

i

wn















(92)

where for the first time we introduced the transpose which, in the matrix convention, swaps columns to rows and rows to columns.
As a general note, a row vector “eats a vector and outputs a scalar”. That is v⊤ : V → R
should be thought of as a 1-form! Due to the matrix conventions, it generally is thecase that “rows eat columns”, that is, every row index should always be thought of as
relating to a 1-form (dual vector), and every column index as relating to a vector. That
is totally consistent to our definition of coordinates.

Maths for Intelligent Systems, Marc Toussaint31

For an orthonormal basis we also havev ∗ (w) = ⟨v, w⟩ = v⊤w .

(93)That is, v⊤ is the coordinate representation of the 1-form v ∗ . (Which also says, thatthe coordinates of the 1-form v ∗ in the special basis (e∗i )ni=1 ⊂ V ∗ coincide with thecoordinates of the vector v.)

3.4

The Structure of Transforms & Singular Value Decomposition

We focus here on linear transforms (or “linear maps”) f : V → U from one vector space
to another (or the same). It turns out that such transforms have a very specific andintuitive structure, which is captured by the singular value decompositon.
3.4.1

The Singular Value Decomposition Theorem

We state the following theorem:
Theorem 3.1 (Singular Value Decomposition). Given two vector spaces V and U
with scalar products, dim(V) = n and dim(U) = m, for every linear transform f :V → U there exist a k ≤ n, m and orthonormal vectors {vi }ki=1 ⊂ V, orthonormal
vectors {ui }ki=1 ⊂ U, and positive scalars σi > 0, i = 1, .., k, such thatf=kX

σi ui vi∗

(94)

i=1

As above, vi∗ = ⟨vi , ·⟩ is the basis 1-form that picks the ith coordinate of a vector
in the basis (vi )ki=1 ⊂ V.aa Note that {v }ki i=1 may not be a full basis of V if k < n. But because {vi } is orthonormal,⟨vi , ·⟩ uniquely picks the ith coordinate no matter how {vi }ki=1 is completed with further n − kvectors to become a full basis.

We first restate this theorem equivalently in coordinates.
Theorem 3.2 (Singular Value Decomposition). For every matrix A ∈ Rm×n there
exists a k ≤ n, m and orthonormal vectors {vi }ki=1 ⊂ Rn , orthonormal vectors{ui } ⊂ Rm , and positive scalars σi > 0, i = 1, .., k, such that
A=

k
X
i=1

⊤
σi ui v⊤i = U SV

(95)

32

Maths for Intelligent Systems, Marc Toussaint

V

U
f

x

span(U )f (x)σspan(V )input (row) space: span(V )input null space: V/span(V )output (column) space: span(U )output null space: U/span(U )PkFigure 2: A linear transformation f = i=1 σi ui v⊤
i can be described as: take the input
x, project it onto the first input fundamental vector v1 to yield a scalar, stretch/squeeze
it by σ1 , and “unproject” this into the first output fundamental vector u1 ; repeat thisfor all i = 1, .., k, and add up the results.
where V = (v1 , .., vk ) ∈ Rn×k , U = (u1 , .., uk ) ∈ Rm×k contain the orthonormal bases vectors as columns and S = diag(σ1 , .., σk ).
Let me rephrase this in a sentence: Every matrix A can be expressed as a linear combination of only k rank-1 matrices. Rank-1 matrices are the most minimalistic kinds ofmatrices and they are always of the form uv⊤ for some u and v. The rank-1 matrix uv⊤
takes an input x, projects it on v (measures its alignment with v), and “unprojects” intou (multiplies v⊤x to the output vector u).
Just to explicitly show the transition from coordinate-free to the coordinate-based theorem, consider arbitrary orthonormal bases {ei }ni=1 ⊂ V and {êi }mi=1 ⊂ U. For x ∈ V wehavef (x) =kX

σi ui ⟨vi , x⟩ =

i=1

=

k
X
i=1

=

k
X

X
X
X
σi (
uli êl )⟨vji ej ,xk ek ⟩

i=1

σi (

X

uli êl )lXvji xk δjk =

jk

i
Xh
U SV ⊤x êl
l

j

l

k
XhX
l

(96)

k

i=1

uli σiXivji xj êl

(97)

j

(98)

l

where vji are the coordinates of vi , uli the coordinates of ui , U = (u1 , .., uk ) is the matrix containing {ui } as columns, V = (v1 , .., vk ) the matrix containing {vi } as columns,
and S = diag(σ1 , .., σk ) the diagonal matrix with elements σi .
We add some definitions based on this:Maths for Intelligent Systems, Marc Toussaint33

Definition 3.14. The rank rank(f ) = rank(A) of a transform f or its matrix A isthe unique minimal k.

Definition 3.15. The determinant of a transform f or its matrix A is
( Qn
± i=1 σi for rank(f ) = n = mdet(f ) = det(A) = det(S) =,0otherwise(99)

where ± depends on whether the transform is a reflection or not.
The last definition is a bit flaky, as the ± is not properly defined. If, alternatively, in theabove theorems we would require V and U to be rotations, that is, elements of SO(n)(of the special
Qn orthogonal group); then negative σ’s would indicate such a reflection anddet(A) = i=1 σi . But above we required σ’s to be strictly positive and V and U only
orthogonal. Fundamental space vectors vi and ui could flip sign. The ± above indicates
how many flip sign.
Definition 3.16. a) The row space (also called right or input fundamental space)rank(f )of a transform f is span{vi }i=1 . The input null space (or right null space) V⊥
is the subspace orthogonal to the row space, such that v ∈ V⊥ ⇒ f (v) = 0.
b) The column space or (also called left or output fundamental space) of a transrank(f )form f is span{ui }i=1 . The output null space (or left null space) U⊥ thesubspace orthogonal to the column space, such that u ∈ U⊥ ⇒ ⟨f (·), u⟩ = 0.

3.5

Point of departure from the coordinate-free notation

The coordinate-free introduction of vectors and transforms helps a lot to understandwhat these fundamentally are. Namely, that coordinate vectors and matrices are ’just’coordinates and rely on a choice of basis; what a metric gij really is; that only for a
Euclidean metric the inner product satisfies ⟨v, w⟩ = v⊤w. Further, the coordinatefree view is essential to understand that vector coordinates behave differently to 1-form
coordinates (e.g., “gradients”!) under a transformation of the basis. We discuss contraversus covariance of gradients at the end of this chapter.
However, we now understood that columns correspond to vectors, rows to 1-forms, andin the Euclidean case the 1-form ⟨v, ·⟩ directly corresponds to v⊤, in the non-Euclidean tov⊤G. In applications we typically represent things from start in orthonormal bases (including perhaps non-Euclidean metrics), there is not much gain sticking to the coordinate-freenotation in most cases. Only when the matrix notation gets confusing (and this happens,e.g. when trying to compute something like the “Jacobian of a Jacobian”, or applying
the chain and product rule for a matrix expression ∂x f (x)⊤A(x)b(x)) it is always a save
harbour to remind what we actually talk about.

34

Maths for Intelligent Systems, Marc ToussaintTherefore, in the rest of the notes we rely on the normal coordinate-based view. Onlyin some explanations we remind at the coordinate-free view when helpful.

3.6

Filling SVD with life

In the following we list some statements—all of them relate to the SVD theorem andtogetherthey’re meant to give a more intuitive understanding of the equation A =
Pk
⊤
⊤
σ
u
i=1 i i vi = U SV .
3.6.1

Understand vv⊤ as a projection

 The projection of a vector x ∈ V onto a vector v ∈ V, is given byx∥ =1v⟨v, x⟩ or
v2

vv⊤
x.
v2

(100)

Here, the v12 is normalizing in case v does not have length |v| = 1.
 The projection-on-v-matrix vv⊤ is symmetric, semi-pos-def, and has rank(vv⊤) =1.
 The projection of a vector x ∈ V onto a subvector space span{vi }ki=1 for orthonormal {vi }ki=1 is given by

x∥ =

X

T
vi v⊤
ix = V V x

(101)

i

where V = (v1 , .., vk ) ∈ Rn×k . The projection matrix V V ⊤ for orthonormal V issymmetric, semi-pos-def, and has rank(V V ⊤) = k.
P The expression i vi v⊤
i is quite related to an SVD. Conversely it shows that the
SVD represent every matrix as a liner combination of kind-of-projects, but thesekind-of-projects uv⊤ first project onto v, but then unproject along u.

3.6.2

SVD for symmetric matrices

 Thm 2 ⇒ Every symmetric matrix A is of the formX⊤A=
λi vi v⊤i = V ΛV(102)

ifor orthonormal V = (v1 , .., vk ). Here λi = ±σi and Λ = diag(λ) is the diagonal
matrix of λ’s. This describes nothing but a stretching/squeezing along orthogonal
projections.

Maths for Intelligent Systems, Marc Toussaint35

 The λi and vi are also the eigenvalues and eigenvectors of A, that is, for alli = 1, .., k:

Avi = λi vi .

(103)

If A has full rank, then the SVD A = V SV ⊤ = V SV -1 is therefore also theeigendecomposition of A.
 The pseudo-inverse of a symmetric matrix is
X
⊤
-1 ⊤
A† =
λ-1
i vi vi = V S V

(104)

i

which simply does the reverse stretching/squeezing along the same orthogonal
projections. Note thatAA† = A† A = V V ⊤(105)

is the projection on {vi }ki=1 . For full rank(A) = n we have V V ⊤ = I andA† = A-1 . For rank(A) < n, we have that A† y minimizes minx ||Ax − y||2 , butthere are infinitly many x’s that minimize this, spanned by the null space of A.
A† y is the minimizer closest to zero (with smallest norm).
n Consider
a data set D = {xi }mi=1 , xi ∈ R . For simplicity assume it has zeroPmmean, i=1 xi = 0. The covariance matrix is defined as

C=

1
1 X
xi x⊤i = X⊤X
n i
n

(106)

where (consistent to ML lecture convention) the data matrix X containes x⊤i in
the ith row. Each xi x⊤i is a projection. C is symmetric and semi-pos-dev. Using
SVD we can writeXC=
λi vi v⊤
(107)
i
i

√
and λi is the data
√ variance along the eigenvector vi ; λi the standard deviation
along vi ; and λi vi the principle axis vectors that make the ellipsoid we typicallyillustrate convariances with.

3.6.3

SVD for general matrices

 For full rank(A) = n, the determinant of a matrix is det(A) = ±define the volume spanned by any {bi }ni=1 as

vol({bi }ni=1 ) = det(B) ,B = (b1 , .., bn ) ∈ Rn×n .

Qi σi . We may

(108)

36

Maths for Intelligent Systems, Marc ToussaintIt follows thatvol({Abi }ni=1 ) = det(A) det(B)(109)

that is, the volume is being multiplied with det(A), which is consistent with ourintuition of transforms as stretchings/squeezings along orthonormal projections.
 The pseudo-inverse of a general matrix isXA† =σi-1 vi u⊤i = V S -1 U⊤ .

(110)

i

If k = n (full input rank), rank(A⊤A) = n and(A⊤A)-1 A⊤ = (V SU⊤U SV ⊤)-1 V SU⊤ = V −⊤S −2 V -1 V SU⊤ = V S -1 U⊤ = A†
(111)
and A† is also called left pseudoinverse because A† A = In .
If k = m (full output rank), rank(AA⊤) = m andA⊤(AA⊤)-1 = V SU⊤(U SV ⊤V SU⊤)-1 = V SU⊤U −⊤S −2 U -1 = V S -1 U⊤ = A†
(112)
and A† is also called right pseudoinverse because AA† = Im .
 Assume m = n (same input/output dimension, or V = U), but k < n. Then there
exist orthogonal V, U ∈ Rn×n such that

A = U DV ⊤ ,S
D = diag(σ1 , .., σk , 0, .., 0) = 

0



0
 .
0

(113)

Here, V and U contain a full orthonormal basis instead of only k orthonormal
vectors. But the diagonal matrix D projects all but k of those to zero. Everysquare matrix A ∈ Rn×n can be written like this.
Definition 3.17 (Rotation). Given a scalar-product ⟨·, ·⟩ on V, a linear transformf : V → V is called rotation iff it preserves the scalar product, that is,∀v, w ∈ V : ⟨f (v), f (w)⟩ = ⟨v, w⟩ .

(114)

 Every rotation matrix is orthogonal, i.e., composed of columns of orthonormal
vectors.
 Every rotation has rank n and σ1,..,n = 1. (No stretching/squeezing.)
 Every square matrix can be written as
rotationU · scalingD · rotationV -1

Maths for Intelligent Systems, Marc Toussaint

3.7

37

Eigendecomposition

Definition 3.18. The eigendecomposition or diagonalization of a square matrix
A ∈ Rn×n is (if it exists!)
A = QΛQ-1(115)

where Λ = diag(λ1 , .., λn ) is a diagonal matrix of eigenvalues. Each column qi of
Q is an eigenvector.
 First note that, unlike SVD, this is not a Theorem but just a definition: If such a
decomposition exists, it is called eigendecomposition. But it exists for almost any
square matrix.
 The set of eigenvalues is the set of roots of the characteristic polynomialpA (λ) = det(λI − A). Why? Because then A − λI has ’volume’ zero (orrank < n), showing that there exists a vector that is mapped to zero, that is
0 = (A − λI)x = Ax − λx.
 Is every square matrix diagonalizable? No, but only an n2 −1-dimensional subset of
matrices are not diagonalizable; most matrices are. A matrix is not diagonalizable
if an eigenvalue λ has multiplicity k (more precisely, λ is a root of pA (λ) withmultiplicity k), but n − rank(A − λI) (the dimensionality of the span of theeigenvectors of λ!) is less than k. Therefore, the eigenvectors of λ are not linearly
independent; they do not span the necessary k dimensions.

So, only very “special” matrices are not diagonalizable. Random matrices are(with prob 1).
 Symmetric matrices? → SVD
 Rotations? Not real. But complex! Think of oscillating projection onto eigenvector. If ϕ is the rotation angle, e±iϕ are eigenvalues.

3.7.1

Power Method

To find the largest eigenvector of A, initialize x randomly and iteratex ← Ax ,

x←

1
x
||x||

(116)

– If this converges, x must be an eigenvector and λ = x⊤Ax the eigenvalue.
– If A is diagonalizable, and x is initially a non-zero linear combination of all eigen
vectors, then it is obvious that x will converge to the “largest” (in absolute terms|λi |) eigenvector (=eigenvector with largest eigenvalue). Actually, if the largest (bynorm) eigenvector is negative, then it doesn’t really converge but flip sign at every
iteration.

38
3.7.2

Maths for Intelligent Systems, Marc ToussaintPower Method including the smallest eigenvalue

A trick, hard to find in the literature, to also compute the smallest eigenvalue and -vector
is the following. We assume all eigenvalues to be positive. Initialize x and y randomly,
iteratex ← Ax ,λ ← ||x|| ,x ← x/λ ,y ← (λI − A)y ,y ← y/||y||(117)

Then y will converge to the smallest eigenvector, and λ − ||y|| will be its eigenvalue.
Note that (in the limit) A − λI only has negative eigenvalues, therefore ||y|| should bepositive. Finding smallest eigenvalues is a common problem in model fitting.
3.7.3

Why should I care about Eigenvalues and Eigenvectors?

– Least squares problems (finding smallest eigenvalue/-vector); (e.g. Camera Calibration, Bundle Adjustment, Plane Fitting)– PCA– stationary distributiuons of Markov Processes– Page Rank
– Spectral Clustering– Spectral Learning (e.g., as approach to training HMMs)

3.8

Beyond this script: Numerics to compute these things

We will not go into details of numerics. Nathan’s script gives a really nice explainationof the QR-method. I just mention two things:(i) The most important forms of matrices for numerics are diagonal matrices, orthogonal
matrices, and upper triangular matrices. One reason is that all three types can very easily
be inverted. A lot of numerics is about finding decompositions of general matrices into
products of these special-form matrices, e.g.:– QR-decomposition: A = QR with Q orthogonal and R upper triangular.
– LU-decomposition: A = LU with U and L⊤ upper triangular.
– Cholesky decomposition: (symmetric) A = C⊤C with C upper triangular– Eigen- & singular value decompositions
Often, these decompositions are intermediate steps to compute eigenvalue or singularvalue decompositions.
(ii) Use linear algebra packages. At the origin of all is LAPACK; browse through http://www.netlib.org/lapack/lug/ to get an impression of what really has been one of
the most important algorithms in all technical areas of the last half century. Modern
wrappers are: Matlab (Octave), which originated as just a console interface to LAPACK;
the C++-library Eigen; or the Python NumPy.

Maths for Intelligent Systems, Marc Toussaint

39

3.9

Derivatives as 1-forms, steepest descent, and the covariant
gradient

3.9.1

The coordinate-free view: A derivative takes a change-of-input vector as
input, and returns a change of output

We previously defined the Jacobian of a function f : Rn → Rd from vector to vector
space. We repeat defining a derivative more generally in coordinate-free form:
Definition 3.19. Given a function f : V → G we define the differential
f (x + hv) − f (x)h→0hdf x : V → G, v 7→ lim(118)

This definition holds whenever G is a continuous space that allows the definition of thislimit and the limit exists (f is differentiable). The notation df |x reads “the differential
at location x”, i.e., evaluating this derivative at location x.
Note that df |x is a mapping from a “tangent vector” v (a change-of-input vector) toan output-change. Further, by this definition df |x is linear. df |x (v) is the directional
derivative we mentioned before. Therefore df |x is a G-valued 1-form. As discussed
earlier, we can introduce coordinates for 1-forms; these coordinates are what typicallyis called the “gradient” or “Jacobian”. But here we explicitly see that we speak of
coordinates of a 1-form.
3.9.2

Contra- and co-variance

In Section 3.2.6 we summarized the effects of a coorindate transformation. We recap
the same here again also for derivatives and scalar products.
We have a vector space V , and a function f : V → R. We’ll be interested in the change
of function value df |x (δ) for change of input δ ∈ V , as well as the value of the scalar
product ⟨δ, δ⟩. All these quantities are defined without any reference to coordinates;
we’ll check now how their coordinate representations change with a change of basis.
As in Section 3.2.6 we have two bases A = (a1 , .., an ) and B = (b1 , .., bn ), and thetransformation T that maps each ai to bi , i.e., B = T A. Given a vector δ we denoteits coordinates in A by [δ]A , and its coordinates in B by [δ]B . Let T AA = B be thematrix representation of T in the old A coordinates (B contains the new basis vectors
b as columns).
 We previously learned that[δ]A = B[δ]B(119)

that is, the matrix B carries new coordinates to old ones. These coordinates are40

Maths for Intelligent Systems, Marc Toussaintsaid to be contra-variant: they transform ‘against’ the transformation of the basis
vectors.
 We require thatdf |x (δ) = [∂x f ]A [δ]A = [∂x f ]B [δ]B

(120)

must be invariant, i.e., the change of function value for δ should not depend onwhether we compute it using A or B coordinates. It follows[∂x f ]A [δ]A = [∂x f ]A B[δ]B = [∂x f ]B [δ]BAB[∂x f ] B = [∂x f ]

(121)
(122)

that is, the matrix B carries old 1-form-coordinates to new 1-form-coordinates.
Therefore, such 1-form-coordinates are called co-variant: they transform ‘with’ thetransformation of basis vectors.
 What we just wrote for the derivative df |x (δ) we could equally write and argue forany 1-form v ∗ ∈ V ∗ ; we always require that the value v ∗ (δ) is invariant.
 We also require that the scalar product ⟨δ, δ⟩ is invariant. Let⟨δ, δ⟩ = [δ]A⊤[G]A [δ]A = [δ]B⊤[G]B [δ]B(123)

where [G]A and [G]B are the 2-form-coordinates (metric tensor) in the old andnew basis. It follows[δ]A⊤[G]A [δ]A = [δ]B⊤B⊤[G]A B[δ]BB⊤

A[G] = B [G] B

(124)
(125)

that is, the matrix carries the old 2-form-coordinates to new ones. These coordinates are called twice co-variant.
Consider the following example: We have the function f : R2 → R, f (x) = x1 +x2 . The
function’s partial derivative is of course ∂f∂x = (1 1). Now let’s transform the coordinatesofthespace:
we
introduce
newcoordinates(z1 , z2 ) = (2x1 , x2 ) or z = B -1 x with B =


1
0
1
2
. The same function, written in the new coordinates, is f (z) =
2 z1 + z2 . The
0 1
1
partial derivatives of that same function, written in these new coordinates, is ∂f∂z = ( 2 1).
Generally, consider we have two kinds of mathematical objects and when we multiply
them together this gives us a scalar. The scalar shouldn’t depend on any choice of
coordinate system and is therefore invariant against coordinate transforms. Then, if
one of the objects transforms in a covariant (“transforming with the transformation”)
manner, the other object must transform in a contra-variant (“transforming contrary
to the transformation”) manner to ensure that the resulting scalar is invariant. This is a
general principle: whenever two things multiply together to give an invariant thing, one
should transform co- the other contra-variant.
Let’s also check Wikipedia:

Maths for Intelligent Systems, Marc Toussaint41– “For a vector to be basis-independent, the components [=coordinates] of the vector
must contra-vary with a change of basis to compensate. That is, the matrix thattransforms the vector of components must be the inverse of the matrix that transforms the basis vectors. The components of vectors (as opposed to those of dual
vectors) are said to be contravariant.
– For a dual vector (also called a covector) to be basis-independent, the componentsof the dual vector must co-vary with a change of basis to remain representing thesame covector. That is, the components must be transformed by the same matrix
as the change of basis matrix. The components of dual vectors (as opposed to thoseof vectors) are said to be covariant.”
Ordinary gradient descent of the form x ← x + α∇f adds objects of different types:
contra-variant coordinates x with co-variant partial derivatives ∇f . Clearly, adding two
such different types leads to an object who’s transformation under coordinate transforms
is strange—and indeed the ordinary gradient descent is not invariant under transformations.
3.9.3

Steepest descent and the covariant gradient vector

Let’s define the steepest descent direction to be the one where, when you make a step oflength 1, you get the largest decrease of f in its linear (=1st order Taylor) approximation.
Definition 3.20. Given f : V → R and a norm ||x||2 = ⟨x, x⟩ (or scalar product)defined on V , we define the steepest descent vector δ ∗ ∈ V as the vector:δ ∗ = argmin df x (δ) s.t. ||δ||2 = 1

(126)

δ

Note that for this definition we need to assume we have a scalar product, otherwise thelength=1 constraint is not defined. Also recall that df |x (δ) = ∂x f (x)δ = ∇f (x)⊤δ areequivalent notations.
Clearly, if we have coordinates in which the norm is Euclidean then||δ||2 = δ⊤δ⇒δ ∗ ∝ −∇f (x)(127)

However, if we have coordinates in which the metric is non-Euclidean, we have:
Theorem 3.3 (Steepest Descent Direction (Covariant gradient)). For a generalscalar product ⟨v, w⟩ = v⊤Gw (with metric tensor G), the steepest descent directionisδ ∗ ∝ −G-1 ∇f (x)(128)

Proof: Let G = B⊤B (Cholesky decomposition) and z = Bδδ ∗ = argmin ∇f⊤δ s.t. δ⊤Gδ = 1
δ

(129)

42

Maths for Intelligent Systems, Marc Toussaint= B -1 argmin ∇f⊤B -1 z s.t. z⊤z = 1

(130)

z
-⊤

∝ B -1 [−B ∇f ] = −G-1 ∇f(131)

For a coordinate transformation B, recall that new metric becomes G̃ = B⊤GB, andff = B⊤∇f . Therefore, the new steepest descent is
the new gradient ∇ff = −B -1 G-1 B -⊤ B⊤∇f = −B -1 G-1 ∇fδe∗ = −[G̃]-1 ∇(132)

and therefore transformes like normal contra-variant coordinates of a vector.
There is an important special case of this, when f is a function over the space of
probability distributions and G is the Fisher metric, which we’ll discuss later.

3.10

Examples and Exercises

3.11

Exercise 1

3.11.1

Basis

Given a linear transform f : R2 → R2 ,
7
f (x) = Ax = 5Consider the basis B =


−10 
x .
−8


 1   2     ,   , which we also simply refer to by the matrix B =
1
1



1 2
2
. Given a vector x in the vector space R , we denote its coordinates in basis B
1 1with xB .
(i) Show that x = BxB .
(ii) What is the matrix F B of f in the basis B, i.e., such that [f (x)]B = F B xB ?
Prove the general equation F B = B -1 AB.
(iii) Provide F B numericallyaNote that for a matrix M = c
3.11.2

−b 
b
1  d
-1
= ad−bc


, M
−c a
d






From the Robotics Course

You have a book lying on the table. The edges of the book define the basis B, the edgesof the table define basis A. Initially A and B are identical (also their origins align). Now
we rotate the book by 45◦ counter-clock-wise about its origin.

Maths for Intelligent Systems, Marc Toussaint43(i) Given a dot p marked on the book at position pB = (1, 1) in the book coordinate
frame, what are the coordinates pA of that dot with respect to the table frame?
(ii) Given a point x with coordinates xA = (0, 1) in table frame, what are its coordinates xB in the book frame?
(iii) What is the coordinate transformation matrix from book frame to table frame,
and from table frame to book frame?
3.11.3

Bases for Polynomials

Consider the set V of all polynomials
for each i = 0, . . . , n.

Pni=0 αi x

iof degree n, where x ∈ R and αi ∈ R(i) Is this set of functions a vector space? Why?
(ii) Consider two different bases of this vector space:
A = {1, x, x2 , . . . , xn }and
B = {1, 1 + x, 1 + x + x2 , . . . , 1 + x + . . . + xn }.
Let f (x) = 1 + x + x2 + x3 be one element in V . (This function f is a vector
in the vector space V , so from here on we refer to it as a vector rather than afunction.)
What are the coordinates [f ]A of this vector in basis A?
What are the coordinates [f ]B of this vector in basis B?
(iii) What matrix I BA allows you to convert between coordinates [f ]A and [f ]B , i.e.
[f ]B = I BA [f ]A ? Which matrix I AB does the same in the opposite direction, i.e.
[f ]A = I AB [f ]B ? What is the relationship between I AB and I BA ?
(iv) What does the difference between coordinates [f ]A − [f ]B represent?
(v) Consider the linear transform t : V → V that maps basis elements of A directly
to basis elements of B:
1→1
x→1+x
x2 → 1 + x + x2..
.
xn → 1 + x + x2 + · · · + xn44

Maths for Intelligent Systems, Marc Toussaint What is the matrix T A for the linear transform t in the basis A, i.e., suchthat [t(f )]A = T A [f ]A ? (Basis A is used for both, input and output spaces.) What is the matrix T B for the linear transform t in the basis B, i.e., suchthat [t(f )]B = T B [f ]B ? (Basis B is used for both, input and output spaces.) What is the matrix T BA if we use A as input space basis, and B as output
space basis, i.e., such that [t(f )]B = T BA [f ]A ?
 What is the matrix T AB if we use B as input space basis, and A as output
space basis, i.e., such that [t(f )]A = T AB [f ]B ?
 Show that T B = I BA T A I AB (cp. Exercise 1(b)). Also note that T AB =T A I AB and T BA = I BA T A .

3.12

Exercise 2

3.12.1

Projections(i) In Rn , a plane (through the origin) is typically described by the linear equationc⊤x = 0 ,(133)

where c ∈ Rn parameterizes the plane. Provide the matrix that describes theorthogonal projection onto this plane. (Tip: Think of the projection as I minus arank-1 matrix.)(ii) In Rn , let’s have k linearly independent {vi }ki=1 , which form the matrix V =(v1 , .., vk ) ∈ Rn×k . Let’s formulate a projection using an optimality principle,
namely,α∗ (x) = argmin ||x − V α||2 .

(134)

α∈Rk

Derive the equation for the optimal α∗ (x) from the optimality principle.
Pk(For information only: Note that V α = i=1 αi vi is just the linear combinationof vi ’s with coefficients α. The projection of a vector x is then x|| = V α∗ (x).)
3.12.2

SVD

Consider the matrices


1 0 0
0 −2 0
A=
0 0 0 ,
0 0 0


1
1
B=
2
0

1
0
1
1


1
−1

0
2

(135)

Maths for Intelligent Systems, Marc Toussaint45(i) Describe their 4 fundamental spaces (dimensionality, possible basis vectors).
(ii) Find the SVD of A (using pen and paper only!)(iii) Given an arbitrary input vector x ∈ R3 , provide the linear transformation matrices
PA and PB that project x into the input null space of matrix A and B, respectively.
(iv) Compute the pseudo inverse A† .
(v) Determine all solutions to the linear equations Ax = y and Bx = y with y =(2, 3, 0, 0). What is the more general expression for an arbitrary y?
3.12.3

Bonus: Scalar product and Orthogonality(i) Show that f (x, y) = 2x1 y1 − x1 y2 − x2 y1 + 5x2 y2 is an scalar product on R2 .
R 2π(ii) In the space of functions with the scalar product ⟨f, g⟩ = 0 f (x)g(x)dx, whatis the scalar product of sin(x) with sin(2x)? (Graphical argument is ok.)(iii) What property does a matrix M have to satisfy in order to be a valid metric tensor,
i.e. such that x⊤M y is a valid scalar product?

3.13

Exercise 3

3.13.1

Eigenvectors(i) A symmetric matrix A ∈ Rn×n is called positive semidefinite (PSD) if x⊤Ax ≥0, ∀x ∈ Rn . (PSD is usually only used with symmetric matrices.) Show that all
eigenvalues of a PSD matrix are non-negative.
(ii) Show that if v is an eigenvector of A with eigenvalue λ, then v is also an eigenvectorof Ak for any positive integer k. What is the corresponding eigenvalue?
(iii) Let v be an eigenvector of A with eigenvalue λ and w an eigenvector of A⊤ witha different eigenvalue µ ̸= λ. Show that v and w are orthogonal.
(iv) Suppose A ∈ Rn×n has eigenvalues λ1 , . . . , λn ∈ R. What are the eigenvalues of
A + αI for α ∈ R and I an identity matrix?
(v) Assume A ∈ Rn×n is diagonalizable, i.e., it has n linearly independent eigenvectors, each with a different eigenvalue. Initialize x ∈ Rn as a random normalized
vector and iterate the two stepsx ← Ax ,

x←

1
x
||x||

46

Maths for Intelligent Systems, Marc ToussaintProve that (under certain conditions) these iterations converge to the eigenvector
x with a largest (in absolute terms |λi |) eigenvalue of A. How fast does thisconverge? In what sense does it converge if the largest eigenvalue is negative?
What if eigenvalues are not different? Other convergence conditions?

(vi) Let A be a positive definite matrix with λmax its largest eigenvalue (in absoluteterms |λi |). What do we get when we apply power iteration method to the matrix
B = A − λmax I? How can we get the smallest eigenvalue of A?
(vii) Consider the following variant of the previous power iteration:z ← Ax ,λ ← x⊤z ,y ← (λI − A)y ,

x←

1
z,
||z||

y←

1
y.
||y||

If A is a positive definite matrix, show that the algorithm can give an estimate of
the smallest eigenvalue of A.
3.13.2

Covariance and PCA

Suppose we’re given a collection of zero-centered data points D = {xi }N
i=1 , with eachxi ∈ Rn . The covariance matrix is defined as
N

C=

1X
1
xi x⊤i = X⊤Xn i=1nwhere (consistent to ML lecture convention) the data matrix X contains each x⊤i as arow, i.e., X⊤ = (x1 , .., xN ).
If we project D onto some unit vector v ∈ Rn , then the variance of the projected data
points is v⊤Cv. Show that the direction that maximizes this variance is the largest
eigenvector of C. (Hint: Expand v in terms of the eigenvector basis of C and exploit
the constraint v⊤v = 1.)
3.13.3

Bonus: RKHS

In machine learning we often work in spaces of functions called Reproducing Kernel
Hilbert Spaces. These spaces are constructed from a certain type of function called thekernel. The kernel k : Rd × Rd → R takes two d-dimensional inputs k(x, x′ ), and fromthe kernel we construct a basis for the space of function, namely B = {k(x, ·)}x∈Rd .
Note that this is a set of infinite element: each x ∈ Rd adds a basis function k(x, ·) to thebasis B. The scalar product between two basis functions kx = k(x, ·) and kx′ = k(x′ , ·)is defined to be the kernel evaluation itself: ⟨kx , kx′ ⟩ = k(x, x′ ). The kernel functionis therefore required to be a positive definite function so that it defines a viable scalar
product.

Maths for Intelligent Systems, Marc Toussaint47(i) Show that for any function f ∈ span B it holds⟨f, kx ⟩ = f (x)(ii) Assume we only have a finite set of points {D = {xi }ni=1 }, which defines a finite
basis {kxi }ni=1 ⊂ B. This finite function basis spans a subspace FD = span{kxi :xi ∈ D} of the space of all functions.
For a general function f , we decompose it f = fs + f⊥ with fs ∈ FD and∀g ∈ FD : ⟨f⊥ , g⟩ = 0, i.e., f⊥ is orthogonal to FD . Show that for every xi ∈ D:
f (xi ) = fs (xi )(Note: This shows that the function values of any function f at the data points
D only depend on the part fs which is inside the spann of {kxi : xi ∈ D}.
This implies the so-called representer theorem, which is fundamental in kernel
machines: A loss can only depend on function values f (xi ) at data points, andtherefore on fs . The part f⊥ can only increase the complexity (norm) of a function.
Therefore, the simplest function to optimize any loss will have f⊥ = 0 and bewithin span{kxi : xi ∈ D}.)(iii) Within span{kxi : xi ∈ D}, what is the coordinate representation of the scalar
product?

4

Optimization

4.1

Downhill algorithms for unconstrained optimization

We discuss here algorithms that have one goal: walk downhill as quickly as possible.
These target at efficiently finding local optima—in constrast to global optimization
methods, which try to find the global optimum.
For such downhill walkers, there are two essential things to discuss: the stepsize andthe step direction. When discussing the stepsize we’ll hit on topics like backtracking
line search, the Wolfe conditions and its implications in a basic convergence proof.
The discussion of the step direction will very much circle around Newton’s method andthereby also cover topics like quasi-Newton methods (BFGS), Gauss-Newton, covariant
and conjugate gradients.
4.1.1

Why you shouldn’t trust the magnitude of the gradientConsider the following 1D function and naive gradient descent x ← x − α∇f for somefixed and small α.

48

Maths for Intelligent Systems, Marc Toussaintsmall gradientsmall step?

large gradientlarge step?

 In plateaus we’d make small steps, at steep slopes (here close to the optimum) we
make huge steps, very likely overstepping the optimum. In fact, for some α thealgorithm might indefinitely loop a non-sensical sequence of very slowly walking
left on the plateau, then accellerating, eventually overstepping the optimum, thenbeing thrown back far to the right again because of the huge negative gradient on
the left.
 Generally, never trust an algorithm that depends on the scaling—or choice ofunits—of a function! An optimization algorithm should be invariant on whether
you measure costs in cents or Euros! Naive gradient descent x ← x − α∇f is not!

As a conclusion, the gradient ∇f gives a reasonable descent direction, but its magnitude
is really arbitrary and no good indication of a good stepsize. Therefore, it often makes
sense to just compute the step directionδ=−1∇f (x)|∇f (x)|(136)

and iterate x ← x + αδ for some appropriate stepsize.

4.1.2

Ensuring monotone and sufficient decrease: Backtracking line search,
Wolfe conditions, & convergence

The first idea is simple: If a step would increase the objective value, reduce the stepsize.
−
We typically use multiplicative stepsize adaptations: Reduce α ← ϱ−
α α with ϱα ≈ 0.5;++
and increase α ← ϱα α with ϱα ≈ 1.2. A simple monotone gradient descent algorithm
reads as follows (the blue part is explained later). Here, the step vector δ is always
normalized and α is adapted on the fly; decreasing when f (x + αδ) is not sufficiently
smaller than f (x).

Maths for Intelligent Systems, Marc Toussaint49

Algorithm 2 Plain gradient descent with backtracking line search
Input: initial x ∈ Rn , functions f (x) and ∇f (x), tolerance θ, parameters (defaults:−ϱ+
α = 1.2, ϱα = 0.5, δmax = ∞, ϱls = 0.01)
1: initialize stepsize α = 1
2: repeat∇f (x)3:δ ← − |∇// (alternative: δ = −∇f (x))
f (x)|⊤
4:
while f (x + αδ) > f (x)+ϱls ∇f (x) (αδ) do// backtracking line search// decrease stepsize5:α ← ϱ−
αα6:
end while
7:x ← x + αδ// increase stepsize
8:
α ← min{ϱ+α α, δmax }9: until |αδ| < θ// perhaps for 10 iterations in sequence

Figure 3: The 1st Wolfe condition: f (x) + ∇f (x)⊤(αδ) is the tangent, which describes
the expected decrease of f (x + αδ) if f was linear. We cannot expect this to be thecase; so f (x + αδ) > f (x) + ϱls ∇f (x)⊤(αδ) weakens this condition.

50

Maths for Intelligent Systems, Marc ToussaintThis sufficiently smaller is described by the blue part and is called the (1st) Wolfe
conditionf (x + αδ) > f (x) + ϱls ∇f (x)⊤(αδ) .

(137)

Figure 3 illustrates this. Note that ∇f (x)⊤(αδ) is a negative value and describes how
much the objective would decrease if f was linear. But f is of course not linear; wecannot expect that a step would really decrease f that much. Instead we require thatit decreases by a fraction of this expectation. ϱls describes this fraction and is typicallychosen very moderate, e.g. ϱls ∈ [0.01, 0.1]. So, the Wolfe conditions requires that f
descreases by the ϱls -fraction of what it would decrease if f was linear. Note that for
α → 0 the Wolfe condition will always be fulfilled for smooth functions, because f“becomes locally linear”.
You’ll proove the following theorem in the exercises. It is fundamental to convex optimization and proves that Alg 2 is efficient for convex objective functions:
Theorem 4.1 (Exponential convergence on convex functions). Let f : Rν → R be
an objective function where the eigenvalues λ of the Hessian ∇2 f (x) are for any x.
Further, assume that Hessian eigenvalues λ are lower bounded by m > 0 and upperbounded by M > m, with m, M ∈ R, at any location x ∈ Rn . (f is convex.) Thenm
Algorithm 2 converges exponentially with convergence rate (1 − 2 Mϱls ϱ−α ) to the
optimum.
But even if your objective function is not globally convex, Alg 2 is an efficient downhill
walker, and once it reaches a convex region it will efficiently walk into its local minimum.
For completeness, there is a second Wolfe condition,|∇f (x + αδ)⊤δ| ≤ b|∇f (x)⊤δ| ,(138)

which states that the gradient magnitude should have decreased sufficiently. We do notuse it much.
4.1.3

The Newton direction

We already discussed the steepest descent direction −G-1 ∇f (x) if G is a metric tensor.
Let’s keep this in mind!
The original Newton method is a method to find the root (that is, zero point) of a′function f (x). In 1D it iterates x ← x − ff′(x)(x) , that is, it uses the gradient f toestimate where the function might cross the x-axis. To find an optimum (minimum
or maximum) of f we want to find the root of its gradient. For x ∈ Rn the Newton
method iteratesx ← x − ∇2 f (x)-1 ∇f (x) .

(139)

Maths for Intelligent Systems, Marc ToussaintNote that the Newton step δ = −∇2 f (x)-1 ∇f (x) is the solution tohi1
min f (x) + ∇f (x)⊤δ + δ⊤∇2 f (x)δ .
δ
2

51

(140)

So the Newton method can also be viewed as 1) compute the 2nd-order Taylor approximation to f at x, and 2) jump to the optimum of this approximation.
Note: If f is just a 2nd-order polynomial, the Newton method will jump to the optimumin just one step.
 Unlike the gradient magnitude |∇f (x)|, the magnitude of the Newton step δ isvery meaningful. It is scale invariant! If you’d rescale f (trade cents by Euros), δis unchanged. |δ| is the distance to the optimum of the 2nd-order Taylor.
 Unlike the gradient ∇f (x), the Newton step δ is truely a vector! The vector itself
is invariant under coordinate transformations; the coordinates of δ transforms
contra-variant, as it is supposed to for vector coordinates.
 The hessian as metric, and the Newton step as steepest descent: Assume
that the hessian H = ∇2 f (x) is pos-def.
Then it fulfils all necessary conditionsP
to define a scalar product ⟨v, w⟩ = ij vi wj Hij , where H plays the role of themetric tensor. If H was the space’s metric, then the steepest descent direction is−H -1 ∇f (x), which is the Newton direction.

Another way to understand the same: In the 2nd-order Taylor approximation f (x+δ) ≈ f (x)+∇f (x)⊤δ+ 21 δ⊤Hδ the Hessian plays the role of a metric tensor. Or: wemay think of the function f as being an isometric parabola f (x + δ) ∝ ⟨δ, δ⟩, but
we’ve chosen coordinates where ⟨v, v⟩ = v⊤Hv and the parabola seems squeezed.
Note that this discussion only holds for pos-dev hessian.
A robust Newton method is the core of many solvers, see Algorithm 3. We do backtracking line search along the Newton direction, but with maximal step size α = 1 (thefull Newton step).
We can additionally add and adapt damping to gain more robustness. Some notes on
the λ:
 In Alg 3, the first line chooses λ to ensure that (∇2 f (x) + λI) is indeed pos-dev—
and a Newton step actually decreases f instead of seeking for a maximum. There
would be other options: instead of adding to all eigenvalues we could only set thenegative ones to some λ > 0.
 δ solves the problemhi11min ∇f (x)⊤δ + δ⊤∇2 f (x)δ + λδ 2 ) .
δ
2
2

(141)

52

Maths for Intelligent Systems, Marc ToussaintAlgorithm 3 Newton method
Input: initial x ∈ Rn , functions f (x), ∇f (x), ∇2 f (x), tolerance θ, parameters (de+−−faults: ϱ+
α = 1.2, ϱα = 0.5, ϱλ = ϱλ = 1, ϱls = 0.01)
1: initialize stepsize α = 1
2: repeat3:choose λ > −(minimal eigenvalue of ∇2 f (x))
4:compute δ to solve (∇2 f (x) + λI) δ = −∇f (x)5:while f (x + αδ) > f (x) + ϱls ∇f (x)⊤(αδ) do// line search// decrease stepsize6:
α ← ϱ−
αα7:
optionally: λ ← ϱ+// increase damping
λ λ and recompute d
8:
end while
9:x ← x + αδ// step is accepted// increase stepsize
10:
α ← min{ϱ+α α, 1}
11:
optionally: λ ← ϱ−// decrease damping
λλ
12: until ||αδ||∞ < θ

So, we added a squared potential λδ 2 to the local 2nd-order Taylor approximation.
This is like introducing a squared penalty for large steps!
 Trust region method: Let’s consider a different mathematical program over thestep:

1
min ∇f (x)⊤δ + δ⊤∇2 f (x)δ s.t. δ 2 ≤ β
δ
2

(142)

This problem wants to find the minimum of the 2nd-order Taylor (like the Newtonstep), but constrained to a stepsize no larger than β. This β defines the trust
region: The region in which we trust the 2nd-order Taylor to be a reasonable
enough approximation.
Let’s solve this using Lagrange parameters (as we will learn it later): Let’s assume
the inequality constraint is active. Then we have1
L(δ, λ) = ∇f (x)⊤δ + δ⊤∇2 f (x)δ + λ(δ 2 − β)2∇δ L(δ, λ) = ∇f (x)⊤ + δ⊤(∇2 f (x) + 2λI)

(143)
(144)

Setting this to zero gives the step δ = −(∇2 f (x) + 2λI)-1 ∇f (x).
Therefore, the λ can be viewed as dual variable of a trust region method. Thereis no analytic relation between β and λ; we cannot determine λ directly fromβ. We could use a constrained optimization method, like primal-dual Newton,
or Augmented Lagrangian approach to solve for λ and δ. Alternatively, we canalways increase λ when the computed steps are too large, and decrease if they aresmaller than β—the Augmented Lagrangian update equations could be used to
do such an update of λ.

Maths for Intelligent Systems, Marc Toussaint53

 The λδ 2 term can be interpreted as penalizing the velocity (stepsizes) of thealgorithm. This is in analogy to a damping, such as “honey pot damping”, inphysical dynamic systems. The parameter λ is therefore also called damping
parameter. Such a damped (or regularized) least-squares method is also called
Levenberg-Marquardt method.
 For λ → ∞ the step direction δ becomes aligned with the plain gradient direction−∇f (x). This shows that for λ → ∞ the Hessian (and metric deformation of thespace) becomes less important and instead we’ll walk orthogonal to the iso-linesof the function.
 The λ term makes the δ non-scale-invariant! δ is not anymore a proper vector!

4.1.4

Least Squares & Gauss-Newton: a very important special caseA special case that appears a lot in intelligent systems is the Least Squares case:
Consider an objective function of the formXf (x) = ϕ(x)⊤ϕ(x) =ϕi (x)2(145)
iwhere we call ϕ(x) the cost features. This is also called a sum-of-squares problem.
We have∂∇f (x) = 2 ∂xϕ(x)⊤ϕ(x)∂∂2∂
ϕ(x)⊤ ∂x
ϕ(x) + 2ϕ(x)⊤ ∂x∇ f (x) = 2 ∂x
2 ϕ(x)
2

(146)
(147)

The Gauss-Newton method is the Newton method for f (x) = ϕ(x)⊤ϕ(x) while approximating ∇2 ϕ(x) ≈ 0. That is, it computes approximate Newton steps∂∂∂
δ = −( ∂xϕ(x)⊤ ∂x
ϕ(x) + λI)-1 ∂xϕ(x)⊤ϕ(x) .

(148)

Note:
∂
∂
 The approximate Hessian 2 ∂xϕ(x)⊤ ∂x
ϕ(x) is always semi-pos-def! Therefore,
no problems arise with negative hessian eigenvalues.

 The approximate Hessian only requires the first-order derivatives of the cost features. There is no need for computationally expensive hessians of ϕ.
 The objective f (x) can be interpreted as the Euclidean norm f (ϕ) = ϕ⊤ϕ but
pulled back into the x-space. More precisely: Consider a mapping ϕ : Rn → Rm
and a general scalar product ⟨·, ·⟩ϕ in the output space. In differential geometry
there is the notion of a pull-back of a metric, that is, we define a scalar product⟨·, ·⟩x in the input space as⟨x, y⟩x = ⟨dϕ(x), dϕ(y)⟩ϕ

(149)

54

Maths for Intelligent Systems, Marc ToussaintAlgorithm 4 Newton method – practical version realized in rai
Input: initial x ∈ Rn , functions f (x), ∇f (x), ∇2 f (x)Input: stopping tolerances θx , θf+−−
Input: parameters ϱ+
α = 1.2, ϱα = 0.5, λ0 = 0.01, ϱλ = ϱλ = 1, ϱls = 0.01
Input: maximal stepsize δmax ∈ R, optionally lower/upper bounds x, x ∈ Rn
initialize stepsize α = 1, λ = λ0
repeat3:
compute smallest eigenvalue σmin of ∇2 f (x) + λI
4:if σmin > 0 is sufficiently posititive then5:
compute δ to solve (∇2 f (x) + λI) δ = −∇f (x)
6:
else
7:
Option 1: λ ← 2λ − σmin and goto line 3// increase regularization
8:
Option 2: δ ← −δmax ∇f (x)/|∇f (x)|// gradient step of length δmax
9:
end if10:
if ||δ||∞ > δmax : δ ← (δmax /||δ||∞ ) δ// cap δ length11:y ← BoundClip(x + αδ, x, x)
12:while f (y) > f (x) + ϱls ∇f (x)⊤(y − x) do// bound-projected line search13:
α ← ϱ−// decrease stepsize
αα14:(unusual option: λ ← ϱ+// increase damping
λ λ and goto line 3)15:y ← BoundClip(x + αδ, x, x)
16:
end while
17:
xold ← x18:x←y// step is accepted19:
α ← min{ϱ+// increase stepsize
α α, 1}20:(unusual option: λ ← ϱ−// decrease damping
λ λ)
21: until ||xold − x||∞ < θx repeatedly, or f (xold ) − f (x) < θf repeatedly
1:
2:

procedure BoundClip(x, x, x)∀i : xi ← min{xi , xi } , xi ← max{xi , xi }24: end procedure
22:

23:

Maths for Intelligent Systems, Marc Toussaint55

where dϕ is the differential of ϕ (a Rm -valued 1-form). Assuming ϕ-coordinates
such that the metric tensor of ⟨·, ·⟩ϕ is Euclidean, we have∂∂⟨x, x⟩x = ⟨dϕ(x), dϕ(x)⟩ϕ = ∂xϕ(x)⊤ ∂x
ϕ(x)

(150)

and therefore, the approximate Hessian is the pullback of a Euclidean costfeature metric, and ⟨x, x⟩x approximates the 2-order polynomial term of f (x),with the non-constant (i.e., Riemannian) pull-back metric ⟨·, ·⟩x .
4.1.5

Quasi-Newton & BFGS: approximating the hessian from gradient observations

To apply full Newton methods we need to be able to compute f (x), ∇f (x), and ∇2 f (x)for any x. However, sometimes, computing ∇2 f (x) is not possible, e.g., because wecannot derive an analytic expression for ∇2 f (x), or it would be too expensive to computethe hessian exactly, or even to store it in memory—especially in very high-dimensional
spaces. In such cases it makes sense to approximate ∇2 f (x) or ∇2 f (x)-1 with a low-rankapproximation.
Assume we have computed ∇f (x1 ) and ∇f (x2 ) at two different points x1 , x2 ∈ Rn . We
definey = ∇f (x2 ) − ∇f (x1 ) ,δ = x2 − x1 .

(151)

From this we may wish to find some approximate Hessian matrix H or H -1 that fulfils
!

H δ=y

or

!

δ = H -1 y

(152)

The first equation is called secant equation. Here are guesses of H and H -1 :
H=

yy⊤
y⊤δ

or

H -1 =

δδ⊤
δ⊤y

(153)

Convince yourself that these choices fulfil the respective desired relation above. However,these choices are under-determined. There exist many alternative H or H -1 that would
be consistent with the observed change in gradient. However, given our understandingof the structure of matrices it is clear that these choices are the lowest rank solutions,
namely rank 1.

Broyden-Fletcher-Goldfarb-Shanno (BFGS): An optimization algorithm computes∇f (xk ) at a series of points x0:K . We incrementally update our guess of H -1 with an
update equationyδ⊤ ⊤ -1 yδ⊤  δδ⊤H I− ⊤ + ⊤ ,H -1 ← I − ⊤
δ y
δ y
δ y

(154)

56

Maths for Intelligent Systems, Marc Toussaintwhich is equivalent to (using the Sherman-Morrison formula)H←H−

Hδδ⊤H⊤ yy⊤+ ⊤ .
δ T Hδ
y δ

(155)

Note:
⊤

 If H -1 is initially zero, this update will assign H -1 ← δδδ⊤y , which is the minimalrank 1 update we discussed above.
-1
 If H -1 is previously non-zero, the red part “deletes certain dimensions” from H .
⊤More precisely, note that I − yδ
y = 0, that is, this rank n − 1 constructionδ⊤ydeletes span{y} from its input space. Therefore, the red part gives zero when
multiplied with y; and it is guaranteed that the resulting H -1 fulfils H -1 y = δ.
The BFGS algorithms uses this H -1 instead of a precise ∇2 f (x)-1 to compute the steps
in a Newton method. All we said about line search and Levenberg-Marquardt damping
is unchanged. 6
In very high-dimensional spaces we do not want to store H -1 densely. Instead we use
aPcompressed storage for low-rank matrices, e.g., storing vectors {vi } such that H -1 =⊤i vi vi . Limited memory BFGS (L-BFGS) makes this more memory efficient: it
limits the rank of the H -1 and thereby the used memory. I do not know the details
myself, but I assume that with every update it might aim to delete the lowest eigenvalue
to keep the rank constant.
4.1.6

Conjugate Gradient

The Conjugate Gradient Method is a method for solving large linear eqn. systems Ax +b = 0. We only mention its extension for optimizing nonlinear functions f (x).
As above, assume that we evaluted ∇f (x1 ) and ∇f (x2 ) at two different points x1 , x2 ∈Rn . But now we make one more assumption: The point x2 is the minimum of a line
search from x1 along the direction δ1 . This latter assumption is quite optimistic: itassumes we did perfect line search. But it gives great information: The iso-lines of f (x)at x2 are tangential to δ1 .
In this setting, convince yourself of the following: Ideally each search direction should beorthogonal to the previous one—but not orthogonal in the conventional Euclidean sense,
but orthogonal w.r.t. the Hessian H. Two vectors d and d′ are called conjugate w.r.t.
a metric H iff d′⊤Hd = 0. Therefore, subsequent search directions should be conjugate
to each other.
6 Taken from Peter Blomgren’s lecture slides:
terminus.sdsu.edu/SDSU/Math693a_f2013/
Lectures/18/lecture.pdf This is the original Davidon-Fletcher-Powell (DFP) method suggested by
W.C. Davidon in 1959. The original paper describing this revolutionary idea – the first quasi-Newton
method – was not accepted for publication. It later appeared in 1991 in the first issue the the SIAM
Journal on Optimization.

Maths for Intelligent Systems, Marc Toussaint57

Conjugate gradient descent does the following:
Algorithm 5 Conjugate gradient descent
Input: initial x ∈ Rn , functions f (x), ∇f (x), tolerance θ
Output: x
1: initialize descent direction d = g = −∇f (x)2: repeat3:
α ← argminα f (x + αd)// line search4:x ← x + αd5:g ′ ← g, gn= −∇f (x) o// store and compute gradβ ← max6:7:8:g⊤(g−g ′ ),0
g ′⊤g ′d ← g + βd
until |∆x| < θ// conjugate descent direction

 The equation for β is by Polak-Ribière: On a quadratic function f (x) = x⊤Hx
this leads to conjugate search directions, d′⊤Hd = 0.
 Intuitively, β > 0 implies that the new descent direction always adds a bit of theold direction. This essentially provides 2nd order information.
 For arbitrary quadratic functions CG converges in n iterations. But this only works
with perfect line search.

4.1.7

Rprop*

Read through Algorithm 6. Notes on this: Stepsize adaptation is done in each coordinate separately !
 The algorithm not only ignores |∇f | but also its exact direction! Only the gradientsigns in each coordinate are relevant. Therefore, the step directions may differ up
to < 90◦ from −∇f .
 It often works surprisingly efficient and robust.
 If you like, have a look at:

Christian Igel, Marc Toussaint, W. Weishui (2005): Rprop using the natural gradient compared to Levenberg-Marquardt optimization. In Trends and Applications
in Constructive Approximation. International Series of Numerical Mathematics,
volume 151, 259-272.

58

Maths for Intelligent Systems, Marc ToussaintAlgorithm 6 Rprop
Input: initial x ∈ Rn , function f (x), ∇f (x), initial stepsize α, tolerance θ
Output: x
1: initialize x = x0 , all αi = α, all gi = 0
2: repeat3:g ← ∇f (x)4:x′ ← x5:
for i = 1 : n do
6:
if gi gi′ > 0 then// same direction as last time7:
αi ← 1.2αi
8:xi ← xi − αi sign(gi )9:gi′ ← gi
10:else if gi gi′ < 0 then// change of direction
11:
αi ← 0.5αi
12:xi ← xi − αi sign(gi )13:
gi′ ← 0// force last case next time
14:
else
15:
xi ← xi − αi sign(gi )16:
gi′ ← gi
17:
end if18:
optionally: cap αi ∈ [αmin xi , αmax xi ]19:
end for′
20: until |x − x| < θ for 10 iterations in sequence

4.2

The general optimization problem – a mathematical program

Definition 4.1. Let x ∈ Rn , f : Rn → R, g : Rn → Rm , h : Rn → Rl . An
optimization problem, or mathematical program, is
min f (x)xs.t. g(x) ≤ 0 ,h(x) = 0

We typically at least assume f, g, h to be differentiable or smooth.
Get an intuition about this problem formulation by considering the following examples.
Always discuss where is the optimum, and at the optimum, how the objective f pulls at
the point, while the constraints g or h push against it.
For the following examples, draw the situation and guess, without much maths, where
the optimum is: A 1D example: x ∈ R, h(x) = sin(x), g(x) = x2 /4 − 1, f some non-linear
function.
 2D example: f (x, y) = x (intuition: the objective is constantly pulling to the left),Maths for Intelligent Systems, Marc Toussaint59

Figure 4: 2D example: f (x, y) = −x, pulling constantly to the right; three inequality
constraints, two active, one inactive. The “pull/push” vectors fulfil the stationarity
condition ∇f + λ1 ∇g1 + λ2 ∇g2 = 0.
h(x, y) = 0 is some non-linear path in the plane → the optimum is at the leftmost tangent-point of h. Tangent-point means that the tangent of h is vertical.
h pushes to the right (always orthogonal to the zero-line of h).
 2D example: f (x, y) = x, g(x, y) = y 2 − x − 1. The zero-line of g is a parabola
towards the right. The objective f pulls into this parabola; the optimum is in the’bottom’ of the parabola, at (−1, 0).
 2D example: f (x, y) = x, g(x, y) = x2 + y 2 − 1. The zero-line of g is a circle.
The objective f pulls to the left; the optimum is at the left tangent-point of thecircle, at (−1, 0).
 Figure 4

4.3

The KKT conditions

Theorem 4.2 (Karush-Kuhn-Tucker conditions). Given a mathematical program,x optimal∇f (x) +⇒∃λ ∈ Rm , κ ∈ Rl s.t.

m
X

l
X

i=1

λi ∇gi (x) +κj ∇hj (x) = 0(stationarity)∀i : gi (x) ≤ 0(primal feasibility)∀i : λi ≥ 0(dual feasibility)j=1∀j : hj (x) = 0 ,∀i : λi gi (x) = 0(complementarity)

Note that these are, in general, only necessary conditions. Only in special cases, e.g.
convex, these are also sufficient.
These conditions should be intuitive in the previous examples:60

Maths for Intelligent Systems, Marc Toussaint The first condition describes the “force balance” of the objective pulling and theactive constraints pushing back. The existance of dual parameters λ, κ could
implicitly be expressed by stating∇f (x) ∈ span({∇g1..m , ∇h1..l })(156)

The specific values of λ and κ tell us, how strongly the constraints push against
the objective, e.g., λi |∇gi | is the force excerted by the ith inequality.
 The fourth condition very elegantly describes the logic of inequality constraints
being either active (λi > 0, gi = 0) or inactive (λi = 0, gi ≤ 0). Intuitively it says:
An inequality can only push at the boundary, where gi = 0, but not inside thefeasible region, where gi < 0. The trick of using the equation λi gi = 0 to express
this logic is beautiful, especially when later we discuss a case which relaxes thisstrict logic to λi gi = −µ for some small µ—which roughly means that inequalities
may push a little also inside the feasible region.
 Special case m = l = 0 (no constraints). The first condition is just the usual∇f (x) = 0.
 Discuss the previous examples as special cases; and how the force balance is met.

4.4

Unconstrained problems to tackle a constrained problem

Assume you’d know about basic unconstrained optimization methods (like standard gradient descent or the Newton method) but nothing about constrained optimization methods. How would you solve a constrained problem? Well, I think you’d very quickly havethe idea to introduce extra cost terms for the violation of constraints—a million people
have had this idea and successfully applied it in practice.
In the following we define a new cost function F (x), which includes the objective f (x)and some extra terms.
Definition 4.2 (Log barrier, squared penalty, Lagrangian, Augmented Lagrangian).
X

Fsp (x; ν, µ) = f (x) + ν

hj (x)2 + µX[gi (x) > 0] gi (x)2j

Flb (x; µ) = f (x) − µ(sqr. penalty)iXlog(−gi (x))(log barrier)i

L(x, λ, κ) = f (x) +Xκj hj (x) +jL̂(x) = f (x) +XjX

λi gi (x)iκj hj (x) +Xi

λi gi (x) +(Lagrangian)Maths for Intelligent Systems, Marc Toussaint61

Figure 5: The function −µ log(−g) (with g on the “x-axis”) for various µ. This isalways undefined (“∞”) for g > 0. For µ → 0 this becomes the hard step function.

+νX

hj (x)2 + µjX[gi (x) > 0] gi (x)2

(Aug. Lag.)i The squared penalty method is straight-forward if we have an algorithm to minimize F (x). We initialize ν = µ = 1, minimize F (x), then increase ν, µ (multiply
with a number > 1) and iterate. For ν, µ → ∞ we retrieve the optimum.
 The log barrier method (see Fig. 5) does exactly the same, except that we decreaseµ towards zero (multiply with a numer < 1 in each iteration). Note that weneed a feasible initialization x0 , because otherwise the barriers are ill-defined! The
whole algorithm will keep the temporary solutions always inside the feasible regions(because the barriers push away from the constraints). That’s why it is also calledinterior point method.
 The Lagrangian is a function L(x, λ, κ) which has the gradient∇L(x, λ, κ) = ∇f (x) +Xi

λi ∇gi (x) +Xκj ∇hj (x) .

(157)

j

That is, ∇L(x, λ, κ) = 0 is our first KKT condition! In that sense, the additionalterms in the Lagrangian generate the push forces of the constraints. If we knew
the correct λ’s and κ’s beforehand, then we could find the optimal x by theunconstrained problem minx L(x, λ, κ) (if this has a unique solution).
 The Augmented Lagrangian L̂ is a function that includes both, squared penalties,
and Lagrangian terms that push proportional to λ, κ. The Augmented Lagrangian
method is an iterative algorithm that, while running, figures out how strongly weneed to push to ensure that the final solution is exactly on the constraints, where
all squared penalties will anyway be zero. It does not need to increase ν, µ andstill converges to the correct solution.

62

Maths for Intelligent Systems, Marc Toussaint4.4.1

Augmented Lagrangian*This is not a main-stream algorithm, but I like it. See (Toussaint, 2014).
In the Augmented Lagrangian L̂, the solver has two types of knobs to tune: the strenghtsof the penalties ν, µ and the strengths of the Lagrangian forces λ, κ. The trick isconceptually easy:
 Initially we set λ, κ = 0 and ν, µ = 1 (or some other constant). In the first
iteration, the unconstrained solver will find x′ = minx L̂(x); the objective f will
typically pull into the penalizations.
 For the second iteration we then choose parameters λ, κ that try to avoid that wewill be pulled into penalizations the next time. Let’s updateκj ← κj + 2µhj (x′ ) ,λi ← max(λi + 2µgi (x′ ), 0).

(158)

Note that 2µhj (x′ ) is the force (gradient) of the equality penalty at x′ ; andmax(λi + 2µgi (x′ ), 0) is the force of the inequality constraint at x′ . What thisupdate does is: it analyzes the forces excerted by the penalties, and translates
them to forces excerted by the Lagrange terms in the next iteration. It tries totrade the penalizations for the Lagrange terms.
More rigorously, observe that, if f, g, h are linear and the same constraints areactive in two consecutive iterations, then this update will guarantee that all penalty
terms are zero in the second iteration, and therefore the solution fulfils the first
KKT condition (Toussaint, 2014). See also the respective exercise.

4.5

The Lagrangian

4.5.1

How the Lagrangian relates to the KKT conditions

The Lagrangian L(x, κ, λ) = f + κ⊤h + λ⊤g has a number of properties that relates itto the KKT conditions:(i) Requiring a zero-x-gradient, ∇x L = 0, implies the 1st KKT condition.
(ii) Requiring a zero-κ-gradient, 0 = ∇κ L = h, implies primal feasibility (the 2nd KKT
condition) w.r.t. the equality constraints.
(iii) Requiring that L is maximized w.r.t. λ ≥ 0 is related to the remaining 2nd and4th KKT conditions:(f (x) if g(x) ≤ 0
max L(x, λ) =
(159)
λ≥0
∞
otherwiseMaths for Intelligent Systems, Marc Toussaint(λ = argmax L(x, λ)

⇒

λ≥0

63

λi = 0if gi (x) < 0(160)
0 = ∇λi L(x, λ) = gi (x) otherwise

This implies either (λi = 0 ∧ gi (x) < 0) or gi (x) = 0, which is equivalent to thecomplementarity and primal feasibility for inequalities.
These three facts show how tightly the Lagrangian is related to the KKT conditions.
To simplify the discussion let us assume only inequality constraints from now on. Fact(i) tells us that if we minx L(x, λ), we reproduce the 1st KKT condition. Fact (iii) tells
us that if we maxλ≥0 L(x, λ), we reproduce the remaining KKT conditions. Therefore,
the optimal primal-dual solution (x∗ , λ∗ ) can be characterized as a saddle point of the
Lagrangian. Finding the saddle point can be written in two ways:
Definition 4.3 (primal and dual problem).
min max L(x, λ)(primal problem)

max min L(x, λ)
x
λ≥0
|
{z
}

(dual problem)

x

λ≥0

l(λ)

(dual function)

Convince
yourself, using 159, thathi the first expression is indeed the original primal problem minx f (x) s.t. g(x) ≤ 0 .

What can we learn from this? The KKT conditions state that, at an optimum, there
exist some λ, κ. This existance statement is not very helpful to actually find them. In
contrast, the Lagrangian tells us directly how the dual parameters can be found: by
maximizing w.r.t. them. This can be exploited in several ways:
4.5.2

Solving mathematical programs analytically, on paper.

Consider the problem
min x2 s.t. x1 + x2 = 1 .

(161)

x∈R2

We can find the solution analytically via the Lagrangian:L(x, κ) = x2 + κ(x1 + x2 − 1)


(162)


1
0 = ∇x L(x, κ) = 2x + κ
 
1

⇒

x1 = x2 = −κ/2

0 = ∇κ L(x, κ) = x1 + x2 − 1 = −κ/2 − κ/2 − 1⇒x1 = x2 = 1/2

⇒

(163)
κ = −1

(164)
(165)

64

Maths for Intelligent Systems, Marc ToussaintHere we first formulated the Lagrangian. In this context, κ is often called Lagrange
multiplier, but I prefer the term dual variable. Then we find a saddle point of L byrequiring 0 = ∇x L(x, κ), 0 = ∇κ L(x, κ). If we want to solve a problem with an inequality
constrained, we do the same calculus for both cases: 1) the constraint is active (handled
like an equality constrained), and 2) the constrained is inactive. Then we check if theinactive case solution is feasible, or the active case is dual-feasible (λ ≥ 0). Note that
if we have m inequality constraints we have to analytically evaluate every combinationof constraints being active/inactive—which are 2m cases. This already hints at the factthat a real difficulty in solving mathematical programs is to find out which inequality
constraints are active or inactive. In fact, if we knew this a priori, everything would
reduce to an equality constrained problem, which is much easier to solve.
4.5.3

Solving the dual problem, instead of the primal.

In some cases the dual function l(λ) = minx L(x, λ) can analytically be derived. In thiscase it makes very much sense to try solving the dual problem instead of the primal.
First, the dual problem maxλ≥0 l(λ) is guaranteed to be convex even if the primal isnon-convex. (The dual function l(λ) is concave, and the constraint λ ≥ 0 convex.)
But note that l(λ) is itself defined as the result of a generally non-convex optimization
problem minx L(x, λ). Second, the inequality constraints of the dual problem are very
simple: just λ ≥ 0. Such inequality constraints are called bound constraints and canbe handled with specialized methods.
However, in general minx maxy f (x, y) ̸= maxy minx f (x, y). For example, in discrete domain x, y ∈ {1, 2}, let f (1, 1) = 1, f (1, 2) = 3, f (2, 1) = 4, f (2, 2) = 2, andminx f (x, y) = (1, 2) and maxy f (x, y) = (3, 4). Therefore, the dual problem is in
general not equivalent to the primal.
The dual function is, for λ ≥ 0, a lower boundhil(λ) = min L(x, λ) ≤ min f (x) s.t. g(x) ≤ 0 .
x

(166)

x

And consequently
max min L(x, λ) ≤ min max L(x, λ)(dual)

x

λ≥0

x

λ≥0

(primal)(167)

We say strong duality holds iff
max min L(x, λ) = min max L(x, λ)
λ≥0

x

x

λ≥0

(168)

If the primal is convex, and there exist an interior point∃x : ∀i : gi (x) < 0(which is called Slater condition), then we have strong duality.

(169)

Maths for Intelligent Systems, Marc Toussaint
4.5.4

65

Finding the “saddle point” directly with a primal-dual Newton method.

In basic unconstrained optimization an efficient way to find an optimum (minimum or
maximum) is to find a point where the gradient is zero with a Newton method. At
saddle points all gradients are also zero. So, to find a saddle point of the Lagrangian we
can equally use a Newton method that seeks for roots of the gradient. Note that sucha Newton method optimizes in the joint primal-dual space of (x, λ, κ).
In the case of inequalities, the zero-gradients view is over-simplified: While facts (i) and(ii) characterize a saddle point in terms of zero gradients, fact (iii) makes this more
precise to handle the inequality case. For this reason it is actually easier to describe theprimal-dual Newton method directly in terms of the KKT conditions: We seek a point(x, λ, κ), with λ ≥ 0, that solves the equation system∇x f (x) + λ⊤∂x g + κ⊤∂x h = 0(170)h(x) = 0(171)

diag(λ)g(x) + µ1m = 0(172)

Note that the first equation is the 1st KKT, the 2nd is the 2nd KKT w.r.t. equalities,
and the third is the approximate 4th KKT with log barrier parameter µ (see below).
These three equations reflect the saddle point properties (facts (i), (ii), and (iii) above).
We define∇f (x) + λ⊤∂g(x) + κ⊤∂h(x) 
r(x, λ, κ) =

h(x)


diag(λ) g(x) + µ1m











(173)

and use the Newton method







x
x
 

 
-1
 
r(x, λ)λ ←  λ  − α ∂r(x, λ, κ)

 
κ
κ






to find the root r(x, λ, κ) = 0 (α is the stepsize). We have

P
P
2
2
2
 ∇ f (x) +i λi ∇ gi (x) +j κj ∇ hj (x)∂r(x, λ, κ) = 

∂h(x)


diag(λ) ∂g(x)(174)∂g(x)⊤∂h(x)⊤


0
0 


diag(g(x))
0
(175)


where ∂r(x, λ, κ) ∈ R(n+m+l)×(n+m+l) . Note that this method uses the hessians∇2 f, ∇2 g and ∇2 h.

The above formulation allows for a duality gap µ. One could choose µ = 0, but oftenthat is not robust. The beauty is that we can adapt µ on the fly, before each Newton
step, so that we do not need a separate outer loop to adapt µ.
Pm1
Before computing a Newton step, we compute the current duality measure µ̃ = − mi=1 λi gi (x1Then we set µ = 2 µ̃ to half of this. In this way, the Newton step will compute a direction

66

Maths for Intelligent Systems, Marc Toussaintthat aims to half the current duality gap. In practise, this leads to good convercence ina single-loop Newton method. (See also Boyd sec 11.7.3.)
The dual feasibility λi ≥ 0 needs to be handled explicitly by the root finder – the line
search can simply clip steps to stay within the bound constraints.
Typically, the method is called “interior primal-dual Newton”, in which case also theprimal feasibility gi ≤ 0 has to be ensured. But I found there are tweaks to make themethod also handle infeasible x, including infeasible initializations.

4.5.5

Log Barriers and the Lagrangian

Finally, let’s revisit the log barrier method. In principle it is very simple: For a given µ,
we use an unconstrained solver to find the minimum x∗ (µ) of
F (x; µ) = f (x) − µXlog(−gi (x)) .

(176)

i

(This process is also called “centering”.) We then gradually decrease µ to zero, always
calling the inner loop to recenter. The generated path of x∗ (µ) is called central path.
The method is simple and has very insightful relations to the KKT conditions and thedual problem. For given µ, the optimality condition is∇F (x; µ) = 0⇒∇f (x) −Xi⇔∇f (x) +Xµ∇gi (x) = 0
gi (x)λi ∇gi (x) = 0 ,(177)
λi gi (x) = −µ

(178)

i

where we defined(!) λi = −µ/gi (x), which guarantees λi ≥ 0 as long as we are in theinterior (gi ≤ 0).
So, ∇F (x; µ) = 0 is equivalent to the modified (=approximate) KKT conditions,where the complemenetarity is relaxed: inequalities may push also inside the feasible
region. For µ → 0 we converge to the exact KKT conditions with strict complementarity.
So µ has the interpretation of a relaxation of complementarity. We can derive anotherinterpretation of µ in terms of suboptimality or the duality gap:
Let x∗ (µ) = minx F (x; µ) be the central path. At each x∗ we may define, as above,
λi = −µ/gi (x∗ ). We note that λ ≥ 0 (dual feasible), as well that x∗ (µ) minimizes the
Lagrangian L(x, λ) w.r.t. x! This is because,
0 = ∇F (x, µ) = ∇f (x) +mX

λi ∇gi (x) = ∇L(x, λ) .

(179)

i=1

Therefore, x∗ is actually the solution to minx L(x, λ), which defines the dual function.

Maths for Intelligent Systems, Marc Toussaint67

We havel(λ) = min L(x, λ) = f (x∗ ) +xmX

λi gi (x∗ ) = f (x∗ ) − mµ .

(180)

i=1

(m is simply the count of inequalities.) That is, mµ is the duality gap between the(suboptimal) f (x∗ ) and l(λ). Further, given that the dual function is a lower bound,l(λ) ≤ p∗ , where p∗ = minx f (x) s.t. g(x) ≤ 0 is the optimal primal value, we havef (x∗ ) − p∗ ≤ mµ .

(181)

This gives the interpretation of µ as an upper bound on the suboptimality of f (x∗ ).

4.6

Convex Problems

We do not put much emphasis on discussing convex problems in this lecture. The
algorithms we discussed so far equally apply on general non-linear programs as well as onconvex problems—of course, only on convex problems we have convergence guarantees,
as we can see from the convergence rate analysis of Wolfe steps based on the assumptionof positive upper and lower bounds on the Hessian’s eigenvalues.
Nevertheless, we at least define standard LPs, QPs, etc. Perhaps the most interesting
part is the discussion of the Simplex algorithm—not because the algorithms is nice orparticularly efficient, but rather because one gains a lot of insights in what actually
makes (inequality) constrained problems hard.
4.6.1

Convex sets, functions, problems

Definition 4.4 (Convex set, convex function). A set X ⊆ V (a subset of some
vector space V ) is convex iff∀x, y ∈ X, a ∈ [0, 1] : ax + (1−a)y ∈ X(182)

A function is defined
convex ⇔ ∀x, y ∈ Rn , a ∈ [0, 1] : f (ax + (1−a)y) ≤ a f (x) + (1−a) f (y)(183)
quasiconvex ⇔ ∀x, y ∈ Rn , a ∈ [0, 1] : f (ax + (1−a)y) ≤ max{f (x), f (y)}(184)
Note: quasiconvex ⇔ for any α ∈ R the sublevel set {x|f (x) ≤ α} is convex.
Further, I call a function unimodal if it has only one local minimum, which is the global
minimum.

68

Maths for Intelligent Systems, Marc ToussaintDefinition 4.5 (Convex program).
Variant 1: A mathematical program minx f (x) s.t. g(x) ≤ 0, h(x) = 0 is convex
if f is convex and the feasible set is convex.
Variant 2: A mathematical program minx f (x) s.t. g(x) ≤ 0, h(x) = 0 is convex
if f and every gi are convex and h is linear.

Variant 2 is the stronger and usual definition. Concerning variant 1, if the feasible set isconvex the zero-level sets of all g’s need to be convex and the zero-level sets of h’s needs
to be linear. Above these zero levels the g’s and h’s could in principle be abribtrarily
non-linear, but these non-linearities are irrelevant for the mathematical program itself.
We could replace such g’s and h’s by convex and linear functions and get the sameproblem.
4.6.2

Linear and quadratic programs

Definition 4.6 (Linear program (LP), Quadratic program (QP)). Special case mathematical programs areLinear Program (LP):min c⊤x s.t. Gx ≤ h, Ax = b

LP in standard form:

min c⊤x s.t. x ≥ 0, Ax = b

Quadratic Program (QP):

x

x

min
x

1 ⊤
x Qx + c⊤x s.t. Gx ≤ h, Ax = b , Qpos-def
2Rarely, also a Quadratically Constrained QP (QCQP) is considered.
An important example for LP are relaxations of integer linear programs,
min c⊤x s.t. Ax = b, xi ∈ {0, 1} ,

(185)

x

which includes Travelling Salesman, MaxSAT or MAP inference problems. Relaxing sucha problem means to instead solve the continuous LP
min c⊤x s.t. Ax = b, xi ∈ [0, 1] .

(186)

x

If one is lucky and the continuous LP problem converges to a fully integer solution, where
all xi ∈ {0, 1}, this is also the solution to the integer problem. Typically, the solutionof the continuous LP will be partially integer (some values converge to the extremexi ∈ {0, 1}, while others are inbetween xi ∈ (0, 1)). This continuous valued solution
gives a lower bound on the integer problem, and provides very efficient heuristics forbacktracking or branch-and-bound search for a fully integer solution.
The standard example for a QP are Support Vector Machines. The primal problem ismin ||β||2 + C
β,ξ

n
X
i=1

ξi s.t. yi (x⊤i β) ≥ 1 − ξi ,ξi ≥ 0(187)

Maths for Intelligent Systems, Marc Toussaint69

the duall(α, µ) = min L(β, ξ, α, µ) = − 14
β,ξ

n X
n
X

αi αi′ yi yi′ x̂⊤i x̂i′ +i=1 i′ =1

n
X

αi

(188)

i=1

max l(α, µ) s.t. 0 ≤ αi ≤ C
α,µ

(189)

(See ML lecture 5:13 for a derivation.)

y

B

A

x

4.6.3

The Simplex Algorithm

Consider an LP. We make the following observations: First, in LPs the equality constraints could be resolved simply by introducing
new coordinates along the zero-hyperplane of h. Therefore, for the conceptual
discussion we neglect equality constraints.
 The objective constantly pulls in the direction −c = −∇f (x).
 If the solution is bounded there need to be some inequality constraints that keep
the solution from travelling to ∞ in the −c direction.
 It follows: The solution will always be located at a vertex, that is, an intersection point of several zero-hyperplanes of inequality constraints.
 In fact, we should think of the feasible region as a polytope that is defined by
all the zero-hyperplanes of the inequalities. The inside the polytope is the feasible
region. The polytope has edges (intersections of two contraint planes), faces, etc.
A solution will always be located at a vertex of the polytope; more precisely, there
could be a whole set of optimal points (on a face orthogonal to c), but at least
one vertex is also optimal.
 An idea for finding the solution is to walk on the edges of the polytope until
an optimal vertex is found. This is the simplex algorithm of Georg Dantzig, 1947.
In practise this procedure is done by “pivoting on the simplex tableaux”—but wefully skip such details here.

70

Maths for Intelligent Systems, Marc Toussaint The simplex algorithm is often efficient, but in worst case it is exponential in both,n and m! This is hard to make intuitive, because the effects of high dimensions are
not intuitive. But roughly, consider that in high dimensions there is a combinatorial
number of ways of how constraints may intersect and form edges and vertices.

Here is a view that much more relates to our discussion of the log barrier method:
Sitting on an edge/face/vertex is equivalent to temporarily deciding which constraints
are active. If we knew which constraints are eventually active, the problem would be
solved: all inequalities become equalities or void. (And linear equalities can directly
be solved for.) So, jumping along vertices of the polytope is equivalent to sequentiallymaking decisions on which constraints might be active. Note though that there are 2m
configurations of active/non-active constraints. The simplex algorithm therefore walks
through this combinatorial space.
Interior point methods do exactly the opposite: Recall that the 4th KKT condition
is λi gi (x) = 0. The log barrier method (for instance) instead relaxes this hard logicof activ/non-active constraints and finds in each iteration a solution to the relaxed 4th
KKT condition λi gi (x) = −µ, which intuitively means that every constraint may be“somewhat active”. In fact, every constraint contributes somewhat to the stationarity
condition via the log barrier’s gradients. Thereby interior point methods
 post-pone the hard decisions about active/non-active constraints
 approach the optimal vertex from the inside of the polytope; avoiding the polytope
surface (and its hard decisions)
 thereby avoids the need to search through a combinatorial space of constraint
activities and instead continuously converges to a decision
 has polynomial worst-case guaranteed

Historically, penalty and barrier methods methods were standard before the Simplex
Algorithm. When SA was discovered in the 50ies, it was quickly considered great. But
then, later in the 70-80ies, a lot more theory was developed for interior point methods,
which now again have become somewhat more popular than the simplex algorithm.

4.6.4

Sequential Quadratic Programming

Just for reference, SQP is another standard approach to solving non-linear mathematicalprograms. In each iteration we compute all coefficients of the 2nd order Taylor f (x+δ) ≈f (x) + ∇f (x)⊤δ + 12 δ⊤Hδ and 1st-order Taylor g(x + δ) ≈ g(x) + ∇g(x)⊤δ and then solve
the QP1
min f (x) + ∇f (x)⊤δ + δ⊤∇2 f (x)δ s.t. g(x) + ∇g(x)⊤δ ≤ 0
δ
2

(190)

Maths for Intelligent Systems, Marc Toussaint71

The optimal δ ∗ of this problem should be seen analogous to the optimal Newton step: If
f were a 2nd-order polynomial and g linear, then δ ∗ would jump directly to the optimum.
However, as this is generally not the case, δ ∗ only gives us a very good direction for line
search. In SQP, we need to backtrack until we found a feasible point and f decreases
sufficiently.
Solving each QP in the sub-routine requires a constrained solver, which itself might have
two nested loops (e.g. using log-barrier or AugLag). In that case, SQP has three nested
loops.

4.7

Blackbox & Global Optimization: It’s all about learning

Even if f, g, h are smooth, the solver might not have access to analytic equations orefficient numeric methods to evaluate the gradients or hessians of these. Therefore we
distinguish (here neglecting the constraint functions g and h):Definition 4.7.
 Blackbox optimization: Only f (x) can be evaluated.
 1st-order/gradient optimization: Only f (x) and ∇f (x) can be evaluated.
 Quasi-Newton optimization: Only f (x) and ∇f (x) can be evaluated, but thesolver does tricks to estimate ∇2 f (x). (So this is a special case of 1st-order
optimization.) Gauss-Newton type optimization: f is of the special form f (x) = ϕ(x)⊤ϕ(x)∂ϕ(x) can be evaluated.
and ∂x
 2nd order optimization: f (x), ∇f (x) and ∇2 f (x) can be evaluated.

In this lecture I very briefly want to add comments on global blackbox optimization.
Global means that we now, for the first time, aim to find the global optimum (withinsome pre-specified bounded range). In essence, to address such a problem we need to
explicitly know what we know about f 7 , and an obvious way to do this is to use Bayesianlearning.

4.7.1

A sequential decision problem formulation

From now on, let’s neglect constraints and focus on the mathematical programmin f (x)x

7 Cf. the KWIK (knows what it knows) framework.

(191)

72

Maths for Intelligent Systems, Marc Toussaintfor a blackbox function f . The optimization process can be viewed as a Markov Decision Process that describes the interaction of the solver (agent) with the function(environment): At step t, Dt = {(xi , yi )}t-1
i=1 is the data that the solver has collected from previous
samples. This Dt is the state of the MDP.
 At step t, the solver may choose a new decision xt about where to sample next.
 Given state Dt and decision xt , the next state is Dt+1 = D ∪ {(xt , f (xt ))}, which
is a deterministic transition given the function f .
 A solver policy is a mapping π : Dt 7→ xt that maps any state (of knowledge) toa new decision.
 We may define an optimal solver policy asZ
π ∗ = argmin⟨yT ⟩ = argmin P (f ) P (DT |π, f ) yT
π

π

(192)

f

where P (DT |π, f ) is deterministic, and P (f ) is a prior over functions.
This objective function cares only about the last value yT sampled by the solver
PT fora fixed time horizon (budget) T . Alternatively, we may choose objectives t=1 yt
PTor t=1 γ t yt for some discounting γ ∈ [0, 1]
The above defined what is an optimal solver! Something we haven’t touched at all
before. The transition dynamics of this MDP is deterministic, given f . However, fromthe perspective of the solver, we do not know f apriori. But we can always computea posterior belief P (f |Dt ) = P (Dt |f ) P (f )/P (Dt ). This posterior belief defines abelief MDP with stochastic transitions
Z Z Z
P (Dt+1 ) =[Dt+1 = D ∪ {(xt , f (xt ))}] π(xt |Dt ) P (f |Dt ) P (Dt ) . (193)
Dt

f

xt

The belief MDP’s state space is P (Dt ) (or equivalently, P (f |Dt ), the current belief overf ). This belief MDP is something that the solver can, in principle, forward simulate—it has all information about it. One can prove that, if the solver could solve its own
belief MDP (find an optimal policy for its belief MDP), then this policy is the optimalsolver policy for the original problem given a prior distribution P (f )! So, in principle we
not only defined what is an optimal solver policy, but can also provide an algorithm tocompute it (Dynamic programming in the belief MDP)! However, this is so expensive tocompute that heuristics need to be used in practise.
One aspect we should learn from this discussion: The solver’s optimal decision is basedon its current belief P (f |Dt ) over the function. This belief is the Bayesian representationof everything one could possibly have learned about f from the data Dt collected so far.
Bayesian Global Optimization methods compute P (f |Dt ) in every step and, based on
this, use a heuristic to choose the next decision.

Maths for Intelligent Systems, Marc Toussaint
4.7.2

73

Acquisition Functions for Bayesian Global Optimization*

In pratice one typically uses a Gaussian Process representation of P (f |Dt ). This meansthat in every iteration we have an estimate fˆ(x) of the function mean and a varianceestimate σ̂(x)2 that describes our uncertainty about the mean estimate. Based on thiswe may define the following acquisition functions
Definition 4.8.
αt (x) =Probability of Improvement (MPI)

R y∗
−∞

N(y|fˆ(x), σ̂(x)) dy(194)

Expected Improvement (EI)
αt (x) =R y∗−∞N(y|fˆ(x), σ̂(x)) (y ∗ − y) dy(195)

Upper Confidence Bound (UCB)
αt (x) = −fˆ(x) + βt σ̂(x)(196)

Predictive Entropy Search Hernández-Lobato et al. (2014)
αt (x) = H[p(x∗ |Dt )] − E{p(y|Dt ; x)} H[p(x∗ |Dt ∪ {(x, y)})](197)= I(x∗ , y|Dt ) = H[p(y|Dt , x)] − E{p(x∗ |Dt )} H[p(y|Dt , x, x∗ )]The last one is special; we’ll discuss it below.
These acquisition functions are heuristics that define how valuable it is to acquire datafrom the site x. The solver then makes the decisionxt = argmax αt (x) .

(198)

x

MPI is hardly being used in practise anymore. EI is classical, originating way back in the
50ies or earlier; Jones et al. (1998) gives an overview. UCB received a lot of attention
recently due to the underlying bandit theory and bounded regret theorems due to thesubmodularity. But I think that in practise EI and UCB perform about equally. As UCB
is somewhat easier to implement and intuitive.
In all cases, note that the solver policy xt = argmaxx αt (x) requires to internally solve
another non-linear optimization problem. However, αt is an analytic function for which
we can compute gradients and hessians which ensures every efficient local convergence.
But again, xt = argmaxx αt (x) needs to be solved globally —otherwise the solver will
also not solve the original problem properly and globally. As a consequence, the optimization of the acquisition function needs to be restarted from many many potential
start points close to potential local minima; typically from grid(!) points over the fulldomain range. The number of grid points is exponential in the problem dimension n.
Therefore, this inner loop can be very expensive.

74

Maths for Intelligent Systems, Marc ToussaintAnd a subjective note: This all sounds great, but be aware that Gaussian Processes withstandard squared-exponential kernels do not generalize much in high dimensions: oneroughly needs exponentially many data points to fully cover the domain and reduce belief
uncertainty globally, almost as if we were sampling from a grid with grid size equal to thekernel width. So, the whole approach is not magic. It just does what is possible given a
belief P (f ). It would be interesting to have much more structured (and heteroscedastic)
beliefs specifically for optimization.
The last acquisition function is called Predictive Entropy Search. This formulation
is beautiful: We sample at places x where the (expected) observed value y informs us
as much as possible about the optimum x∗ of the function! Formally, this means tomaximize the mutual information between y and x∗ , in expectation over y|x.
4.7.3

Classical model-based blackbox optimization (non-global)*

A last method very worth mentioning: Classical model-based blackbox optimization
simply fits a local polynomial model to the recent data and takes this a basis for search.
This is similar to BFGS, but now for the blackbox case where we not even observe
gradients. See Algorithm 7.
The local fitting of a polynomial model is again a Machine Learning method. Whether
this gives a function approximation for optimization depends on the quality of the data Dtused for this approximation. Classical model-based optimization has interesting heuristics to evaluate the data quality as well as sample new points to improve the data
quality. Here is a rough algorithm (following Nodecal et al.’s section on “Derivative-free
optimization”):
Some notes:
 Line 4 implement an explicit trust region approach, which hard bound α on thestep size.
 Line 5 is like the Wolfe condition. But here, the expected decrease is [f (x̂) −fˆ(x̂ + δ)] instead of −αδ∇f (x).
 If there is no sufficient decrease we may blame it on two reasons: bad data or a
too large stepsize.
 Line 10 uses the data determinant as a measure of quality! This is meant in thesense of linear regression on polynomial features. Note that, with data matrix
X ∈ Rn×dim(β) , β̂ ls = (X⊤X)-1 X⊤y is the optimal regression. The determinantdet(X⊤X) or det(X) = det(D) is a measure for well the data supports theregression. If the determinant is zero, the regression problem is ill-defined. The
larger the determinant, the lower the variance of the regression estimator.
 Line 11 is an explicit exploration approach: We add a data point solely for the
purpose of increasing the data determinant (increasing the data spread). Interest-

Maths for Intelligent Systems, Marc Toussaint75

Algorithm 7 Classical model-based optimization
1
1: Initialize D with at least
2 (n + 1)(n + 2) data points
2: repeat
3:Compute a regression fˆ(x) = ϕ2 (x)⊤β on D
4:Compute δ = argminδ fˆ(x̂ + δ) s.t. |δ| < α
5:
if f (x̂ + δ) < f (x̂) − ϱls [f (x̂) − fˆ(x̂ + δ)] then// test sufficient decrease
6:
Increase the stepsize α
7:Accept x̂ ← x̂ + δ
8:
Add to data, D ← D ∪ {(x̂, f (x̂))}9:else// no sufficient decrease10:
if det(D) is too small then// blame the data quality
11:
Compute x+ = argmaxx′ det(D ∪ {x′ }) s.t. |x − x′ | < α
12:
Add to data, D ← D ∪ {(x+ , f (x+ ))}13:
else// blame the stepsize
14:
Decrease the stepsize α
15:
end if16:
end if17:
Perhaps prune the data, e.g., remove argmaxx∈∆ det(D \ {x})18: until x converges

ing. Nocedal describes in more detail a geometry-improving procedure to update
D.

4.7.4

Evolutionary Algorithms*

There are interesting and theoretically well-grounded evolutionary algorithms for optimization, such as Estimation-of-Distribution Algorithms (EDAs). But generally, don’t
use them as first choice.

4.8

Examples and Exercises

4.9

Exercise 6

4.9.1

Convergence proofa) Given a function f : Rn → R with fmin = minx f (x). Assume that its Hessian—that is, the eigenvalues of ∇2 f —are lower bounded by m > 0 and upper bounded by
M > m, with m, M ∈ R. Prove that for any x ∈ Rn it holds
f (x) −

1
1|∇f (x)|2 ≤ fmin ≤ f (x) −|∇f (x)|2 .
2m
2M

76

Maths for Intelligent Systems, Marc ToussaintTip: Start with bounding f (x) between the functions with maximal and minimal curvature. Then consider the minima of these bounds. Note, it also follows:|∇f (x)|2 ≥ 2m(f (x) − fmin ) .
b) Consider backtracking line search with Wolfe parameter ϱls ≤ 12 , and step decreaseϱ−1αfactor ϱ−
α . First prove that line search terminates the latest when M ≤ α ≤ M , andthen it found a new point y for whichf (y) ≤ f (x) −ϱls ϱ−α|∇f (x)|2 .
M

From this, using the result from a), prove the convergence equationhi2mϱls ϱ−αf (y) − fmin ≤ 1 −(f (x) − fmin ) .
M
4.9.2

Backtracking Line SearchConsider the functions
fsq (x) = x⊤Cx ,(199)

fhole (x) = 1 − exp(−x⊤Cx) .

(200)
i−1

with diagonal matrix C and entries C(i, i) = c n−1 , where n is the dimensionality of x.
We choose a conditioning8 c = 10. To plot the function for n = 2, you can use gnuplot
calling

set isosamples 50,50
set contourf(x,y) = x*x+10*y*y#f(x,y) = 1 - exp(-x*x-10*y*y)
splot [-1:1][-1:1] f(x,y)a) Implement gradient descent with backtracking, as described on page 42 (Algorithm
2 Plain gradient descent). Test the algorithm on fsq (x) and fhole (x) with start pointx0 = (1, 1). To judge the performance, create the following plots:– The function value over the number of function evaluations.
– For n = 2, the function surface including algorithm’s search trajectory. If using
gnuplot, store every evaluated point x and function value f (x) in a line (with n + 1
entries) in a file ’path.dat’, and plot using
unset contour
splot [-3:3][-3:3] f(x,y), ’path.dat’ with lines

8 The word “conditioning” generally denotes the ratio of the largest and smallest Eigenvalue of the
Hessian.

Maths for Intelligent Systems, Marc Toussaint77b) Play around with parameters. How does the performance change for higher dimensions, e.g., n = 100? How does the performance change with ρls (the Wolfe stop
criterion)? How does the alternative in step 3 work?
c) Newton step: Modify the algorithm simply by multiplying C -1 to the step. How doesthat work?
(The Newton direction diverges (is undefined) in the concave part of fhole (x). We’recheating here when always multiplying with C -1 to get a good direction.)

4.9.3

Gauss-Newton

In x ∈ R2 consider the function
f (x) = ϕ(x)⊤ϕ(x) ,sin(ax1 )sin(acx2 )

ϕ(x) =  2x1 2cx2

The function is plotted above for a = 4 (left) and a = 5 (right, having local minima),
and conditioning c = 1. The function is non-convex.
Extend your backtracking method implemented in the last week’s exercise to a GaussNewton method (with constant λ) to solve the unconstrained minimization problem
minx f (x) for a random start point in x ∈ [−1, 1]2 . Compare the algorithm for a = 4
and a = 5 and conditioning c = 3 with gradient descent.
4.9.4

Robust unconstrained optimization

A ’flattened’ variant of the Rosenbrock function is defined as
f (x) = log[1 + (x2 − x21 )2 +1
(1 − x2 )2 ]100

and has the minimum at x∗ = (1, 1). For reference, the gradient and hessian are
g(x) := 1 + (x2 − x21 )2 +1
(1 − x2 )2
100

(201)

78

Maths for Intelligent Systems, Marc Toussainti
1 h− 4(x2 − x21 )x1
g(x)
i
1 h
2
∂x2 f (x) =2(x2 − x21 ) −(1 − x2 )
g(x)
100
i
h
i2
1 h 2
8x1 − 4(x2 − x21 )∂x21 f (x) = − ∂x1 f (x) +
g(x)
h
i2
1 h
2 i
∂x22 f (x) = − ∂x2 f (x) +
2+
g(x)
100
h
ih
i
i
1 h
∂x1 ∂x2 f (x) = − ∂x1 f (x) ∂x2 f (x) +
− 4x1
g(x)
∂x1 f (x) =

(202)
(203)
(204)
(205)
(206)

a) Use gnuplot to display the function copy-and-pasting the following lines:

set isosamples 50,50
set contourf(x,y) = log(1+(y-(x**2))**2 + .01*(1-x)**2 ) - 0.01
splot [-3:3][-3:4] f(x,y)(The ’-0.01’ ensures that you can see the contour at the optimum.) List and discuss atleast three properties of the function (at different locations) that may raise problems tonaive optimizers.
b) Use x = (−3, 3) as starting point for an optimization algorithm. Try to code an
optimization method that uses all ideas mentioned in the lecture. Try to tune it to beefficient on this problem (without cheating, e.g. by choosing a perfect initial stepsize.)

4.10

Exercise 7

4.10.1

Lagrangian Method of Multipliersc
In a previous exercise we defined the “hole function” fhole(x). Assume conditioningc = 10 and use the Lagrangian Method of Multipliers to solve on paper the following
constrained optimization problem in 2D:c
min fhole(x) s.t. h(x) = 0xh(x) = v⊤x − 1

(207)
(208)

Near the very end, you won’t be able to proceed until you have special values for v. Go
as far as you can without the need for these values.

Maths for Intelligent Systems, Marc Toussaint4.10.279

Equality Constraint Penalties and Augmented Lagrangian

The squared penalty approach to solving a constrained optimization problem minimizes
min f (x) + µxmXhi (x)2 .

(209)

i=1

The Augmented Lagrangian method adds a Lagrangian term and minimizes
min f (x) + µxmXhi (x)2 +

i=1

m
X

λi hi (x) .

(210)

i=1

Assume that we first minimize (209) we end up at a minimum xb.
Now prove that setting λi = 2µhi (bx) will, if we assume that the gradients ∇f (x) and∇h(x) are (locally) constant, ensure that the minimum of (210) fulfills the constraints
h(x) = 0.

4.11

Exercise 8

4.11.1

Lagrangian and dual function(Taken roughly from ‘Convex Optimization’, Ex. 5.1)
Consider the optimization problem
min x2 + 1xs.t.

(x − 2)(x − 4) ≤ 0

with variable x ∈ R.
a) Derive the optimal solution x∗ and the optimal value p∗ = f (x∗ ) by hand.
b) Write down the Lagrangian L(x, λ). Plot (using gnuplot or so) L(x, λ) over x forvarious values of λ ≥ 0. Verify the lower bound property minx L(x, λ) ≤ p∗ , where p∗
is the optimum value of the primal problem.
c) Derive the dual function l(λ) = minx L(x, λ) and plot it (for λ ≥ 0). Derive the dual
optimal solution λ∗ = argmaxλ l(λ). Is maxλ l(λ) = p∗ (strong duality)?
4.11.2

Optimize a constrained problem

Consider the following constrained problem
min
x

n
X
i=1

xi s.t. g(x) ≤ 0

(211)

80

Maths for Intelligent Systems, Marc Toussaint⊤

x x − 1g(x) = 


−x1

(212)

a) First, assume x ∈ R2 is 2-dimensional, and draw on paper what the problem looks
like and where you expect the optimum.
b) Find the optimum analytically using the Lagrangian. Here, assume that you know
apriori that all constraints are active! What are the dual parameters λ = (λ1 , λ2 )?
Note: Assuming that you know a priori which constraints are active is a huge assumption!
In real problems, this is the actual hard (and combinatorial) problem. More on this laterin the lecture.
c) Implement a simple the Log Barrier Method. Tips:– Initialize x = ( 21 , 12 ) and µ = 1– First code an inner loop:– In each iteration, first compute the gradient of the log-barrier function. RecallthatX
F (x; µ) = f (x) − µlog(−gi (x))
(213)
i

X
∇F (x; µ) = ∇f − µ(1/gi (x))∇gi (x)

(214)

i

– Then perform a backtracking line search along −∇F (x, µ). In particular, backtrack if a step goes beyond the barrier (where g(x) ̸≤ 0 and F (x, µ) = ∞).
– Iterate until convergence; let’s call the result x∗ (µ). Further, compute λ∗ (m) =−(µ/g1 (x), µ/g2 (x)) at convergence.
– Decrease µ ← µ/2, recompute x∗ (µ) (with the previous x∗ as initialization) anditerate this.
Does x∗ and λ∗ converge to the expected solution?
Note: The path x∗ (µ) = argminx F (x; µ) (the optimum in dependence of µ) is called
central path.

4.12

Exercise 9

Comment: Solving problems in the real world involves 2 parts:
1) formulating the problem as an optimization problem (conform to a standard optimization problem category) (→ human)2) the actual optimization problem (→ algorithm)

Maths for Intelligent Systems, Marc Toussaint81

These exercises focus on the first type, which is just as important as the second, as itenables the use of a wider range of solvers. Exercises from Boyd et al http://www.
stanford.edu/~boyd/cvxbook/bv_cvxbook.pdf:
4.12.1

Network flow problem

Solve Exercise 4.12 (pdf page 193) from Boyd & Vandenberghe, Convex Optimization.
4.12.2

Minimum fuel optimal control

Solve Exercise 4.16 (pdf page 194) from Boyd & Vandenberghe, Convex Optimization.
4.12.3

Reformulating an ℓ1 -norm(This is a subset of Exercise 4.11 (pdf page 193) from Boyd & Vandenberghe.)n
Let x ∈ RP
. The optimization problem is minx ||x − b||1 , where the ℓ1 -norm is definednas ||z||1 = i=1 |zi |. Reformulate this optimization problem as a Linear Program.

4.13

Exercise 10

4.13.1

Restarts of Local Optima

The following function is essentially the Rastrigin function, but written slightly differently.
It can be tuned to become uni-modal and is a sum-of-squares problem. For x ∈ R2 wedefinesin(ax1 )sin(acx2 )

f (x) = ϕ(x)⊤ϕ(x) , ϕ(x) = 
 2x1 2cx2
The function is plotted above for a = 4 (left) and a = 5 (right, having local minima),
and conditioning c = 1. The function is non-convex.

82

Maths for Intelligent Systems, Marc ToussaintChoose a = 6 or larger and implement a random restart method: Repeat initializingx ∼ U([−2, 2]2 ) uniformlly, followed by a gradient descent (with backtracking line search
and monotone convergence).
Restart the method at least 100 times. Count how often the method converges to which
local optimum.

4.13.2

GP-UCB Bayesian Optimization

Find an implementation of Gaussian Processes for your language of choice (e.g. python:scikit-learn, or Sheffield/Gpy; octave/matlab: gpml) and implement GP-UCB global
optimization. Test your implementation with different hyperparameters (Find the best
combination of kernel and its parameters in the GP) on the 2D function defined above.
On the webpage you find a starting code to use GP regression in scikit-learn. To installscikit-learn: https://scikit-learn.org/stable/install.html

5Probabilities & Information

It is beyond the scope of these notes to give a detailed introduction to probability theory.
There are excellent books:
 Thomas & Cover Bishop MacKay

Instead, we first recap very basics of probability theory, that I assume the reader has
already seen before. The next section will cover this. Then we focus on specific topics that, in my opinion, deepen the understanding of the basics, such as the relation
between optimization and probabilities, log-probabilities & energies, maxEntropy andmaxLikelihood, minimal description length and learning.

5.1

BasicsFirst, in case you wonder about justifications of the use of (Bayesian) probabilities versus
fuzzy sets or alike, here some pointers to look up: 1) Cox’s theorem, which derives frombasic assumptiond about “rationality and consistency” the standard probability axioms;
2) t-norms, which generalize probability and fuzzy calculus; and 3) read about objective
vs. subjecte Bayesian probability.

Maths for Intelligent Systems, Marc Toussaint
5.1.1

83

Axioms, definitions, Bayes ruleDefinition 5.1 (set-theoretic axioms of probabilities).
 An experiment can have multiple outcomes; we call the set of possible outcomes sample space or domain S
 A mapping P : A ⊆ S 7→ [0, 1], that maps any subset A ⊆ S to a real
number, is called probability measure on S iff– P (A) ≥ 0 for any A ⊆ S (non-negativity)SP– P ( i Ai ) = i P (Ai ) if Ai ∩ Aj = ∅ (additivity)– P (S) = 1 (normalization)
 Implications are:– 0 ≤ P (A) ≤ 1– P (∅) = 0– A ⊆ B ⇒ P (A) ≤ P (B)– P (A ∪ B) = P (A) + P (B) − P (A ∩ B)
– P (S \ A) = 1 − P (A)Formally, a random variable X is a mapping X : S → Ω from a measureable space S(that is, a sample space S that has a probability measure P ) to another sample spaceΩ, which I typically call the domain dom(X) of the random variable. Thereby, themapping X : S → Ω now also defines a probability measure over the domain Ω:
P (B ⊆ Ω) = P ({s : X(s) ∈ B})(215)

In practise we just use the following notations:Definition 5.2 (Random Variable).
 Let X be a random variable with discrete domain dom(X) = Ω P (X = x) ∈ R denotes the specific probability that X = x for some x ∈ Ω
 P (X) denotes the probability distribution (function over Ω) ∀x∈Ω : 0 ≤ P (X = x) ≤ 1
P
x∈Ω P (X = x) = 1

PP We often use the short hand X P (X) · · · = x∈dom(X) P (X = x) · · ·when summing over possible values of a RV
If we have two or more random variables, we have

84

Maths for Intelligent Systems, Marc ToussaintDefinition 5.3 (Joint, marginal, conditional, independence, Bayes’ Theorem).
 We denote the joint distribution of two RVs as P (X, Y )P The marginal is defined as P (X) = Y P (X, Y )) The conditional is defined as P (X|Y ) = PP(X,Y(Y ) , which fulfilsP∀Y :X P (X|Y ) = 1.

 X is independent of Y iff P (X, Y ) = P (X) P (Y ), or equivalently,
P (X|Y ) = P (X).
 The definition of a conditional implies the product rule

P (X, Y ) = P (X|Y ) P (Y ) = P (Y |X) P (X)(216)

and Bayes’ Theorem
P (X|Y ) =P (Y |X) P (X)P (Y )(217)

The individual terms in Bayes’ Theorem are typically given names:
posterior =likelihood · priornormalization(218)(Sometimes, the normalization is also called evidence.)
 X is conditionally independent of Y given Z iff P (X|Y, Z) = P (X|Z) or
P (X, Y |Z) = P (X|Z) P (Y |Z)

5.1.2

Standard discrete distributions

BernoulliRVx ∈ {0, 1}

parameterµ ∈ [0, 1]

Betaµ ∈ [0, 1]α, β ∈ R+

Multinomialx ∈ {1, .., K}

Dirichletµ ∈ [0, 1]K , ||µ||1 =1µ ∈ [0, 1]K , ||µ||1 =1
α1 , .., αK ∈ R+

distribution
Bern(x | µ) = µx (1 −µ)1−x
Beta(µ | a, b)
=
1
a−1
b−1
µ
(1
−
µ)
B(a,b)P (x = k | µ) = µk
Dir(µ | α) ∝

QK

αk −1
k=1 µkClearly, the Multinomial is a generalization of the Bernoulli, as the Dirichlet is of theBeta. The mean of the Dirichlet is ⟨µi ⟩ = Pαiαj , its mode is µ∗i = Pααi −1. The modej −K
jof a distribution p(x) is defined as argmaxx p(x).

jMaths for Intelligent Systems, Marc Toussaint
5.1.3

85

Conjugate distributions

Definition 5.4 (Conjugacy). Let p(D|x) be a likelihood conditional on a RV x. A
family C of distributions (i.e., C is a space of distributions, like the space of all Beta
distributions) is called conjugate to the likelihood function p(D|x) iff
p(x) ∈ C⇒

p(x|D) =p(D|x) p(x)
∈C.
p(D)

The standard conjugates you should know:
RVlikelihoodµBinomial Bin(D | µ)µMultinomial Mult(D | µ)µGauss N(x | µ, Σ)
λ
1D Gauss N(x | µ, λ-1 )Λ
nD Gauss N(x | µ, Λ-1 )(µ, Λ)nD Gauss N(x | µ, Λ-1 )

5.1.4

(219)

conjugate
Beta Beta(µ | a, b)
Dirichlet Dir(µ | α)Gauss N(µ | µ0 , A)
Gamma Gam(λ | a, b)
Wishart Wish(Λ | W, ν)
Gauss-Wishart
N(µ | µ0 , (βΛ)-1 ) Wish(Λ | W, ν)

Distributions over continuous domain

Definition 5.5. Let x be a continuous RV. The probability density function (pdf)p(x) ∈ [0, ∞) defines the probability
Z b
P (a ≤ x ≤ b) =p(x) dx ∈ [0, 1]

(220)

a

The cumulative probability distribution F (y) = P (x ≤ y) =[0, 1] is the cumulative integral with limy→∞ F (y) = 1

Ry−∞

dx p(x) ∈However, I and most others say probability distribution to refer to probability
density function.
One comment about integrals. If p(x) is a probability density function and f (x) somearbitrary function, typically one writes
Zf (x) p(x) dx ,(221)xwhere dx denotes the (Borel) measure we integrate over. However, some authors (correctly) think of a distribution p(x) as being a measure over the space dom(x) (insteadof just a function). So the above notation is actually “double” w.r.t. the measures. So
they might (also correctly) writeZp(x) f (x) ,(222)
x

and take care that there is exactly one measure to the right of the integral.

86

Maths for Intelligent Systems, Marc Toussaint5.1.5

Gaussian

Definition 5.6. We define an n-dim Gaussian in normal form as
N(x | µ, Σ) =1
1
exp{− (x − µ)⊤ Σ-1 (x − µ)}1/22| 2πΣ |

(223)

with mean µ and covariance matrix Σ. In canonical form we define
N[x | a, A] =exp{− 21 a⊤A-1 a}1exp{− x⊤ A x + x⊤a}-1
1/2
2
| 2πA |

(224)

with precision matrix A = Σ-1 and coefficient a = Σ-1 µ (and mean µ = A-1 a).
Gaussians are used all over—below we explain in what sense they are the probabilistic analogue to a parabola (or a 2nd-order Taylor expansions). The most important
properties are:
 Symmetry:N(x | a, A) = N(a | x, A) = N(x − a | 0, A)

 Product:
N(x | a, A) N(x | b, B) = N[x | A-1 a + B -1 b, A-1 + B -1 ] N(a | b, A + B)
N[x | a, A] N[x | b, B] = N[x | a + b, A + B] N(A-1 a | B -1 b, A-1 + B -1 )
 “Propagation”:
R
N(x | a + F y, A) N(y | b, B) dy = N(x | a + F b, A + F BF⊤)y
 Transformation:N(F x + f | a, A) = | F1 | N(x | F -1 (a − f ), F -1 AF -⊤ )
 Marginal& conditional:


A C
x a
N
,
= N(x | a, A) · N(y | b + C⊤A-1 (x - a), B − C⊤A-1 C)
y b C⊤ B

More Gaussian identities are found at http://ipvs.informatik.uni-stuttgart.de/mlr/marc/
notes/gaussians.pdf

Example 5.1 (ML estimator of the mean of a Gaussian). Assume we have data
D = {x1 , .., xn }, each xi ∈ Rn , with likelihood
Q
P (D | µ, Σ) = i N(xi | µ, Σ)
(225)
n
1X
argmax P (D | µ, Σ) =
xi
(226)
n i=1
µ
n

argmax P (D | µ, Σ) =Σ1X(xi − µ)(xi − µ)⊤
n i=1

(227)

Maths for Intelligent Systems, Marc Toussaint87

Assume we are initially uncertain about µ (but know Σ). We can express thisuncertainty using again a Gaussian N[µ | a, A]. Given data we haveQ
P (µ | D) ∝ P (D | µ, Σ) P (µ) = i N(xi | µ, Σ) N[µ | a, A]
(228)
Q
P
= i N[µ | Σ-1 xi , Σ-1 ] N[µ | a, A] ∝ N[µ | Σ-1 i xi , nΣ-1 + A]
(229)
Note: in the limit A → 0 (uninformative prior) this becomes
P (µ | D) = N(µ |11X
xi , Σ)
n i
n

(230)

which is consistent with the Maximum Likelihood estimator

5.1.6“Particle distribution”

Usually, “particles” are not listed as standard continuous distribution. However I thinkthey should be. They’re heavily used in several contexts, especially as approximatingother distributions in Monte Carlo methods and particle filters.
Definition 5.7 (Dirac or δ-distribution). In distribution theory it is proper to definea distribution δ(x) that is the derivative of the Heavyside step function H(x),∂δ(x) = ∂x
H(x) ,

H(x) = [x ≥ 0] .

(231)

It is akward to think of δ(x) as a normal function, as it’d be “infinite” at zero. But
at least we understand that is has the propertiesZδ(x) = 0 everywhere except at x = 0 ,δ(x) dx = 1 .
(232)

I sometimes call the Dirac distribution also a point particle: it has all its unit “mass”
concentrated at zero.
Definition 5.8 (Particle Distribution). We define a particle distribution q(x) asa mixture of Diracs,
q(x) :=NX

wi δ(x − xi ) ,

(233)

i=1
n
which is parameterized by the number N , the locations {xi }N
i=1 , xi ∈ R , and theNnormalized weights {wi }i=1 , wi ∈ R, ||w||1 = 1 of the N particles.

We say that a particle distribution q(x) approximates another distribution p(x) iff for88

Maths for Intelligent Systems, Marc Toussaintany (smooth) fZ⟨f (x)⟩p =f (x)p(x)dx ≈x

PN

i=1 wi f (xi )(234)Note the generality of this statement! f could be anything, it could be any features ofthe variable x, like coordinates of x, or squares, or anything. So basically this statement
says, whatever you might like to estimate about p, you can approximate it based on theparticles q.
Computing particle approximations of comples (non-analytical, non-tracktable) distributions p is a core challenge in many fields. The true p could for instance be a distributions
over games (action sequences). The approximation q could for instance be samples generated with Monte Carlo Tree Search (MCTS). The tutorial An Introduction to MCMC
for Machine Learning www.cs.ubc.ca/~nando/papers/mlintro.pdf gives an excellent introduction. Here are some illustrations of what it means to approximate somep by particles q, taken from this tutorial. The black line is p, histograms illustrate theparticles q by showing how many of (uniformly weighted) particles fall into a bin:(from de Freitas et al.)5.2

Between probabilities and optimization: neg-log-probabilities,
exp-neg-energies, exponential family, Gibbs and Boltzmann

There is a natural relation between probabilities and “energy” (or “error”). Namely, ifp(x) denotes a probability for every possible value of x, and E(x) denotes an energy forstate x—or an error one assigns to choosing x—then a natural relation isp(x) = e−E(x) ,E(x) = − log p(x) .

(235)

Maths for Intelligent Systems, Marc Toussaint89

Why is that? First, outside the context of physics it is perfectly fair to just define
axiomatically an energy E(x) as neg-log-probability. But let me try to give some morearguments for why this is a useful definition.
Let assume we have p(x). We want to find a quantity, let’s call it error E(x), which
is a function of p(x). Intuitively, if a certain value x1 is more likely than another,p(x1 ) > p(x2 ), then picking x1 should imply less error, E(x1 ) < E(x2 ) (Axiom 1).
Further, when we have two independent random variables x and y, probabilities are
multiplicative, p(x, y) = p(x)p(y). We require axiomatically that error is additive,E(x, y) = E(x) + E(y). From both follows that E needs to be some logarithm of p!
The same argument, now more talking about energy : Assume we have two independent(physical) systems x and y. p(x, y) = p(x)p(y) is the probability to find them in certain
states. We axiomatically require that energy is additive, E(x, y) = E(x) + E(y).
Again, E needs to be some logarithm of p. In the context of physics, what could
be questioned is “why is p(x) a function of E(x) in the first place?”. Well, that ismuch harder to explain and really is a question about statistical physics. Wikipedia
under keywords “Maxwell-Boltzmann statistics” and “Derivation from microcanonicalensemble” gives an answer. Essentially the argument is a follows: Given many many
molecules
Pn in a gas, each of which can have a different energy ei . The total energy
E = i=1 ei must be conserved. What is the distribution over energy levels that has
the most microstates? The answer is the Boltzmann distribution. (And why do we,in nature, find energy distributions that have the most microstates? Because these aremost likely.)
Bottom line is: p(x) = e−E(x) , probabilities are multiplicative, energies or errors additive.
Let me state some fact just to underline how useful this way of thinking is: Given an energy function E(x), its Boltzmann distribution is defined asp(x) = e−E(x) .

(236)This is sometimes also called Gibbs distribution.
 In machine learning, when data D is given and we have some model β, we typicallytry to maximize the likelihood p(D|β). This is equivalent to minimizing the neglog-likelihood

L(β) = − log p(D|β) .

(237)

This neg-log-likelihood is a typical measure for error of the model. And this error
is additive w.r.t. the data, whereas the likelihood is multiplicative, fitting perfectly
to the above discussion.
 The Gaussian distribution p(x) ∝ exp{− 12 ||x − µ||2 /σ 2 } is related to the error
E(x) = 12 ||x − µ||2 /σ 2 , which is nothing but the squared error with the precision
matrix as metric. That’s why squared error measures (classical regression) and
Gaussian distributions (e.g., Bayesian Ridge regression) are directly related.

A Gaussian is the probabilistic analoque to a parabola.

90

Maths for Intelligent Systems, Marc Toussaint The exponential family is defined asp(x|β) = h(x)g(β) exp{β⊤ϕ(x)}(238)

Often h(x) = 1, so let’s neglect this for now. The key point is that the energy
is linear in the features ϕ(x). This is exactly how discriminative functions (forclassification in Machine learning) are typically formulated.
In the continuous case, the features ϕ(x) are often chosen as basis polynomials—
just as in polynomial regression. Then, β are the coefficients of the energy polynomial and the exponential family is just the probabilistic analogue to the spaceof polynomials.
 When we have many variables x1 , .., xn , the structure of a cost function overthese
P variables can often be expressed as being additive in terms: f (x1 , .., xn ) =The respectivei ϕi (x∂i ) where ∂i denotes the ith group of variables.
Q
PBoltzmann
distribution is a factor graph p(x1 , .., xn ) ∝ i fi (x∂i ) = exp{ i βi ϕi (x∂i )where ∂i denotes theSo, factor graphs are the probabilistic analoque to additive functions.
 − log p(x) is also the “optimal” coding length you should assign to a symbol x.
PEntropy is expected error: H[p] = x −p(x) log p(x) = ⟨− log p(x)⟩p(x) , wherep itself it used to take the expectation.

Assume you use a “wrong” distribution q(x) to decide on the Rcoding length of symbols drawn from p(x). The expected length of an encoding is x p(x)[− log q(x)] ≥H(p).
The Kullback-Leibler divergence is the difference:Zp(x)D p q =
p(x) log≥0q(x)xProof of inequality, using the Jenson inequality:
Z
Z
q(x)
q(x)
− p(x) log≥ − log p(x)
=0
p(x)
p(x)
x
x

(239)

(240)

So, my message is that probabilities and error measures are naturally related. However,in the first case we typically do inference, in the second we optimize. Let’s discuss therelation between inference and optimization a bit more. For instance, given data D andparameters β, we may define
Definition 5.9 (ML, MAP, and Bayes estimate). Given data D and a parametericmodel p(D|β), we define Maximum likelihood (ML) parameter estimate:β ML := argmaxβ P (D|β)

Maths for Intelligent Systems, Marc Toussaint91

 Maximum a posteriori (MAP) parameter estimate:

β MAP = argmaxβ P (β|D)
 Bayesian parameter estimate:

P (β|D) ∝ P (D|β) P (β)
used for Bayesian prediction: P (prediction|D) =RβP (prediction|β) P (β|D)

Both, the MAP and the ML estimates are really just optimization problems.
The Bayesian parameter estimate P (β|D), which can then be used to do fully Bayesianprediction, is in principle different. However, in practise also here optimization is acore tool for estimating such distributions if they cannot be given analytically. This is
described next.

5.3

Information, Entropie & Kullback-Leibler

Consider the following problem. We have data drawn i.i.d. from p(x) where x ∈ X insome discrete space X. Let’s call every x a word. The problem is to find a mappingfrom words to codes, e.g. binary codes c : X → {0, 1}∗ . The optimal solution is inprinciple simple: Sort all possible words in a list, ordered by p(x) with more likely words
going first; write all possible binary codes in another list, with increasing code lengths.
Match the two lists, and this is the optimal encoding.
Let’s try to get a more analytical grip of this: Let l(x) = |c(x)| be the actual code
length assigned to word x, which is an integer value. Let’s defineq(x) =

1 −l(x)
2
Z

(241)

P
with the normalization constrant Z = x 2−l(x) . Then we haveXX
X
p(x)[− log2 q(x)] = −p(x) log 2−l(x) +p(x) log Z
x

x∈X

=

X

p(x)l(x) + log Z .

(242)

x

(243)

x

What about log Z? Let l-1 (s) be the set of words that have been assigned codes oflength l. There can only be a limited number of words encoded with a given length. For
instance, |L-1 (1)| must not be greater than 2, |L-1 (2)| must not be greater than 4, and|l-1 (s)| must not be greater than 2l . We haveX∀s :[l(x) = s] ≤ 2s
(244)
x∈X

∀s :

X
x∈X

[l(x) = s]2−s ≤ 1

(245)

92

Maths for Intelligent Systems, Marc Toussaint∀s :X

2−l(x) ≤ 1

(246)

x∈X

However, this way of thinking is ok for separated codes. If such codes would be in acontinuous stream of bits you’d never know where a code starts or ends. Prefix codes
fix this problem by defining a code tree with leaves that clearly define when a code ends.
For prefix codes it similarly holds
XZ=
2−l(x) ≤ 1 ,
(247)
x∈X

which is called Kraft’s inequality. That finally givesXX
p(x)[− log2 q(x)] ≤p(x)l(x)

5.4

(248)

x

x∈X

The Laplace approximation: A 2nd-order Taylor of log p

Assume we want to estimate some q(x) we cannot express analytically. E.g., q(x) =p(x|D) ∝ P (D|x)p(x) for some awkward likelihood function p(D|x). An example fromrobotics is: x is stochastically controlled path of a robot. p(x) is a prior distribution overpaths that includes how the robot can actually move and some Gaussian prior (squared
costs!) over controlls. If the robot is “linear”, p(x) can be expressed nicely and analytically; it if it non-linear, expressing p(x) is already hard. However, p(D|x) might indicatethat we do not see collisions on the path—but collisions are a horrible function, usually
computed by some black-box collision detection packages that computes distances between convex meshes, perhaps giving gradients but certainly not some analytic function.
So q(x) can clearly not be expressed analytically.
One way to approximate q(x) is the Laplace approximation
Definition 5.10 (Laplace approximation). Given a smooth distribution q(x), we
define its Laplace approximation asq̃(x) = exp{−Ẽ(x)} ,(249)

where Ẽ(x) is the 2nd-order Taylor expansion
1Ẽ(x) = E(x∗ ) + (x − x∗ )⊤∇2 E(x∗ )(x − x∗ )
2

(250)

of the energy E(x) = − log q(x) at the modex∗ = argmin E(x) = argmax q(x) .
x

(251)

x

First, we observe that the Laplace approximation is a Gaussian, because its energy isa parabola. Further, notice that in the Taylor expansion we skipped the linear term.
That’s because we are at the mode x∗ where ∇E(x∗ ) = 0.

Maths for Intelligent Systems, Marc Toussaint93

The Laplace approximation is the probabilistic analoque of a local second-order approximation of a function, just as we used it in Newton methods. However, it is defined tobe taken specifically at the mode of the distribution.
Now, computing x∗ is a classical optimization problem x∗ = argminx E(x) which one
might ideally solve using Newton methods. These Newton methods anyway computethe local Hessian of E(x) in every step—at the optimum we therefore have the Hessianalready, which is then the precision matrix of our Gaussian.
The Laplace approximation is nice, very efficient to use, e.g., in the context of optimal
control and robotics. While we can use the expressive power of probability theory toformalize the problem, the Laplace approximation brings us computationally back toefficient optimization methods.

5.5

Variational Inference

Another reduction of inference to optimization is variational inference.

Definition 5.11 (variational inference). Given a distribution p(x), and a parameterized family of distributions q(x|β), the variational approximation of p(x) is definedas
argmin D q p
(252)
q

5.6

The Fisher information metric: 2nd-order Taylor of the KLD

Recall our notion of steepest descent—it depends on the metric in the space!
Consider the space of probability distributions p(x; β) with parameters β. We think ofevery p(x; β) as a point in the space and wonder what metric is useful to compare twopoints p(x; β1 ) and p(x; β2 ). Let’s take the KLD
TODO: Let p ∈ ΛX , that is, p is a probability distribution over the space X. Further, letθ ∈ Rn and θ 7→ p(θ) is some parameterization of the probability distribution. Thenthe derivative dθ p(θ) ∈ Tp ΛX is a vector in the tangent space of ΛX . Now, for suchvectors, for tangent vectors of the space of probability distributions, there is a generic
metric, the Fisher metric: [TODO: move to ’probabilities’ section]

94

Maths for Intelligent Systems, Marc Toussaint5.7

Examples and Exercises

5.8

Exercise 10

Note: These exercises are for ’extra credits’. We’ll discuss them on Thu, 21th Jan.

5.8.1

Maximum Entropy and Maximum Likelihood(These are taken from MacKay’s book Information Theory..., Exercise 22.12 & .13)a) Assume that a random variable x with discrete domain dom(x) = X comes from aprobability distribution of the form
P (x | w) =

d
i
hX
1
wk fk (x) ,
exp
Z(w)
k=1

where the functions fk (x) are given, and the parameters w ∈ Rd are not known. A
data set D = {xi }ni=1 ofPn points x is supplied. Show by differentiating the log
nlikelihood log P (D|w) =i=1 log P (xi |w) that the maximum-likelihood parametersw∗ = argmaxw log P (D|w) satisfynX

P (x | w∗ ) fk (x) =

x∈X

1 X
fk (xi )n i=1

where the left-hand sum is over all x, and the right-hand sum is over the data points.
A shorthand for this result is that each function-average under the fitted model must
equal the function-average found in the data:⟨fk ⟩P (x | w∗ ) = ⟨fk ⟩Db) When confronted by a probability distribution P (x) about which only a few facts are
known, the maximum entropy principle (MaxEnt) offers a rule for choosing a distributionthat satisfies those constraints. According to MaxEnt, you should select the P (x) thatmaximizes the entropy
X
H(P ) = −P (x) log P (x)xsubject to the constraints. Assuming the constraints assert that the averages of certain
functions fk (x) are known, i.e.,⟨fk ⟩P (x) = Fk ,

Maths for Intelligent Systems, Marc Toussaint95

show, by introducing Lagrange multipliers (one for each constraint, including normalization), that the maximum-entropy distribution has the form
hXi
1
PMaxEnt (x) =exp
wk fk (x)Zkwhere the parameters Z and wk are set such that the constraints are satisfied. And hence
the maximum entropy method gives identical results to maximum likelihood fitting ofan exponential-family model.

5.9

Exercise 11

Note: The exercise will take place on Tue, 2nd Feb. Hung will also prepare how much
‘votes’ you collected in the exercises.

5.9.1

Maximum likelihood and KL-divergence

Assume we have a very large data set D = {xi }ni=1 of samples xi ∼ q(x) from somedata distribution q(x). Using this data set we can approximate any expectation⟨f ⟩q =

X

q(x)f (x) ≈xnXf (xi ) .

i=1

Assume we have a parameteric family of distributions p(x|β) and would find the Maximum Likelihood (ML) parameter β ∗ = argmaxβ p(D|β). Express this ML problem as a
KL-divergence minimization.
5.9.2

Laplace Approximation

In the context of so-called “Gaussian Process Classification” the following problem arises(we neglect dependence on x here): We have a real-valued RV f ∈ R with prior P (f ) =N(f | µ, σ 2 ). Further we have a Boolean RV y ∈ {0, 1} with conditional probability
P (y = 1 | f ) = σ(f ) =ef.
1 + ef

The function σ is called sigmoid funtion, and f is a discriminative value which predicts
y = 1 if it is very positive, and y = 0 if it is very negative. The sigmoid function has theproperty∂σ(f ) = σ(f ) (1 − σ(f )) .
∂f

96

Maths for Intelligent Systems, Marc ToussaintGiven that we observed y = 1 we want to compute the posterior P (f | y = 1), whichcannot be expressed analytically. Provide the Laplace approximation of this posterior.

(Bonus)RAs an alternative to the sigmoid function σ(f ), we can use the probit functionzϕ(z) = −∞ N(x|0, 1) dx to define the likelihood P (y = 1 | f ) = ϕ(f ). Now how can
the posterior P (f | y = 1) be approximated?
5.9.3

Learning = Compression

In a very abstract sense, learning means to model the distribution p(x) for given data
D = {xi }ni=1 . This is literally the case for unsupervised learning; regression, classificationand graphical model learning could be viewed as specific instances of this where x factoresin several random variables, like input and output.
Show in which sense the problem of learning is equivalent to the problem of compression.
5.9.4

A gzip experiment

Get three text files from the Web, approximately equal length, mostly text (no equations
or stuff). Two of them should be in English, the third in Frensh. (Alternatively, perhaps,
not sure if it’d work, two of them on a very similar topic, the third on a very different.)
How can you use gzip (or some other compression tool) to estimate the mutual information between every pair of files? How can you ensure some “normalized” measures
which do not depend too much on the absolute lengths of the text? Do it and check
whether in fact you find that two texts are similar while the third is different.

(Extra) Lempel-Ziv algorithms (like gzip) need to build a codebook on the fly. How doesthat fit into the picture?

5.9.5

Maximum Entropy and ML(These are taken from MacKay’s book Information Theory..., Exercise 22.12 & .13)a) Assume that a random variable x with discrete domain dom(x) = X comes from aprobability distribution of the form
P (x | w) =

d
hX
i
1
exp
wk fk (x) ,
Z(w)

(253)

k=1

where the functions fk (x) are given, and the parameters w ∈ Rd are not known. A
data set D = {xi }ni=1 of n points x is supplied. Show by differentiating the log

Maths for Intelligent Systems, Marc Toussaint97

Pnlikelihood log P (D|w) =i=1 log P (xi |w) that the maximum-likelihood parametersw∗ = argmaxw log P (D|w) satisfynX

P (x | w∗ ) fk (x) =

x∈X

1 X
fk (xi )
n i=1

(254)

where the left-hand sum is over all x, and the right-hand sum is over the data points.
A shorthand for this result is that each function-average under the fitted model must
equal the function-average found in the data:⟨fk ⟩P (x | w∗ ) = ⟨fk ⟩D(255)b) When confronted by a probability distribution P (x) about which only a few facts are
known, the maximum entropy principle (MaxEnt) offers a rule for choosing a distributionthat satisfies those constraints. According to MaxEnt, you should select the P (x) thatmaximizes the entropy
X
H(P ) = −P (x) log P (x)(256)
xsubject to the constraints. Assuming the constraints assert that the averages of certain
functions fk (x) are known, i.e.,⟨fk ⟩P (x) = Fk ,(257)

show, by introducing Lagrange multipliers (one for each constraint, including normalization), that the maximum-entropy distribution has the form
hXi
1
PMaxEnt (x) =exp
wk fk (x)
(258)
Z
k

where the parameters Z and wk are set such that the constraints are satisfied. And hence
the maximum entropy method gives identical results to maximum likelihood fitting ofan exponential-family model.

A

Gaussian identities

Definitions
We define a Gaussian over x with mean a and covariance matrix A as the function
N(x | a, A) =11
exp{− (x-a)⊤ A-1 (x-a)}
1/2
2
|2πA|

(259)

with property N (x | a, A) = N (a| x, A). We also define the canonical form with precision
matrix A as
N[x | a, A] =exp{− 12 a⊤A-1 a}1exp{− x⊤ A x + x⊤a}2|2πA-1 |1/2

(260)

98

Maths for Intelligent Systems, Marc Toussaintwith properties
N[x | a, A] = N(x | A-1 a, A-1 )(261)

N(x | a, A) = N[x | A a, A ] .

(262)

-1

-1

Non-normalized GaussianN(x, a, A) = |2πA|1/2 N(x|a, A)
1= exp{− (x-a)⊤ A-1 (x-a)}
2

(263)
(264)

Matrices [matrix cookbook: http://www.imm.dtu.dk/pubdb/views/edoc_download.php/3274/
pdf/imm3274.pdf](A-1 + B -1 )-1 = A (A+B)-1 B = B (A+B)-1 A
-1

-1 -1

-1

(A − B ) = A (B-A)

(265)

B

(266)

∂x |Ax | = |Ax | tr(A-1x ∂x Ax )-1-1∂x Ax = −Ax (∂x Ax ) A-1
x
-1
-1
-1

(267)
(268)
-1

-1

-1

(A + U BV ) = A − A U (B + V A U ) V A
-1

-1 -1

-1

-1

(A + B ) = A − A(B + A) A
⊤

-1 ⊤

⊤

-1

-1 ⊤

-1

(270)
-1 ⊤ -1(A + J BJ) J B = A J (B + JA J )
⊤

(269)

-1 ⊤

(A + J BJ) A = I − (A + J BJ) J BJ

(271)
(272)

(269)=Woodbury; (271,272) holds for pos def A and B
Derivatives∂x N(x|a, A) = N(x|a, A) (−h⊤) ,h := A-1 (x-a)∂θ N(x|a, A) = N(x|a, A) ·hi
11
− h⊤(∂θ x) + h⊤(∂θ a) − tr(A-1 ∂θ A) + h⊤(∂θ A)h2
2
h 1
1 ⊤ -1⊤∂θ N[x|a, A] = N[x|a, A] − x ∂θ Ax + a A ∂θ AA-1 a
2
2
i
1
⊤
⊤ -1
-1
+ x ∂θ a − a A ∂θ a + tr(∂θ AA )2
∂θ Nx (a, A) = Nx (a, A) ·hi
1h⊤(∂θ x) + h⊤(∂θ a) + h⊤(∂θ A)h2
Product
The product of two Gaussians can be expressed as
N(x | a, A) N(x | b, B)

(273)

(274)

(275)

(276)

Maths for Intelligent Systems, Marc Toussaint99= N[x | A-1 a + B -1 b, A-1 + B -1 ] N(a | b, A + B) ,

(277)

-1

= N(x | B(A+B) a + A(A+B) b, A(A+B) B) N(a | b, A + B) ,(278)

= N[x | a + b, A + B] N(A-1 a | B -1 b, A-1 + B -1 )

(279)

-1

-1

N[x | a, A] N[x | b, B]= N[x| . . . ] N[A a | A(A+B) b, A(A+B) B]
-1

-1

-1

(280)

= N[x| . . . ] N[A a | (1-B(A+B) ) b, (1-B(A+B) ) B] ,
-1

-1

-1

(281)

N(x | a, A) N[x | b, B]= N[x | A-1 a + b, A-1 + B] N(a | B -1 b, A + B -1 )(282)

= N[x| . . . ] N[a | (1-B(A +B) ) b, (1-B(A +B) ) B]
-1

-1

-1

-1

Convolution
R
N(x | a, A) N(y − x | b, B) dx = N(y | a + b, A + B)
x

(283)

(284)

DivisionN(x|a, A)N(x|b, B) = N(x|c, C)

-1

-1



N(c|b, C + B)

-1

C c=A a−B b
C -1 = A-1 − B -1
N[x|a, A] N[x|b, B] ∝ N[x|a − b, A − B]

(285)
(286)

Expectations
Let x ∼ N(x | a, A),
E{x} g(x) :=RxN(x | a, A) g(x) dx⊤

E{x} x = a ,

(287)

⊤

E{x} xx = A + aa(288)

E{x} f + F x = f + F a
⊤

(289)

⊤

E{x} x x = a a + tr(A)
⊤

(290)
⊤

E{x} (x-m) R(x-m) = (a-m) R(a-m) + tr(RA)Transformation(291)

Linear transformations imply the following identities,

N(x | a, A) = N(x + f | a + f, A) ,N(x | a, A) = |F | N(F x | F a, F AF⊤)
(292)

N(F x + f | a, A) ==1
N(x | F -1 (a − f ), F -1 AF -⊤ )|F |1
N[x | F⊤A-1 (a − f ), F⊤A-1 F ] ,
|F |

(293)
(294)

100

Maths for Intelligent Systems, Marc ToussaintN[F x + f | a, A] =“Propagation”
respectively)

1
N[x | F⊤(a − Af ), F⊤AF ] .
|F |

(295)

(propagating a message along a coupling, using eqs (277) and (283),yN(x | a + F y, A) N(y | b, B) dy = N(x | a + F b, A + F BF⊤)yN(x | a + F y, A) N[y | b, B] dy = N[x | (F

R
R

-⊤

-1

(296)

-K)(b + BF a), (F-⊤K = F -⊤ B(F -⊤ A-1 F -1 +B)-1-K)BF -1 ] ,
(297)
(298)

marginal & conditional:


a
A
x
,
y b + Fa FA

A⊤F⊤
B +F A⊤F⊤N(x | a, A) N(y | b + F x, B) = N
(299)


A C
x a
,
N
= N(x | a, A) · N(y | b + C⊤A-1 (x-a), B − C⊤A-1 C)
y b C⊤ B(300)

x a + F⊤B -1 b A + F⊤B -1 F −F⊤B -1
N[x | a, A] N(y | b + F x, B) = N,y−B -1 F
B -1
B -1 b
(301)


⊤ -1
⊤ -1
⊤
x a + F B b A + F B F −F
N[x | a, A] N[y | b + F x, B] = N
,
y
b
−F
B
(302)


A C
x a
N
,
= N[x | a − CB -1 b, A − CB -1 C⊤] · N[y | b − C⊤x, B]
y b C⊤ B
(303)
-1
b
A C
b = |A|b |B| , where A = A − CB D
= |A| |B|b
D B
B = B − DA-1 C
(304)


A
D

C
B

"

-1
=

"
=

b-1
A
-1
b
−B DA-1

b -1−A-1 C B
-1
b
B

#

b-1
A
b-1
−B -1 DAb-1 CB -1
−A
b -1
B

#

(305)

(306)

pair-wise belief We have a message α(x) = N[x|s, S], transition P (y|x) = N(y|Ax +a, Q), and a message β(y) = N[y|v, V ], what is the belief b(y, x) = α(x)P (y|x)β(y)?
b(y, x) = N[x|s, S] N(y|Ax + a, Q-1 ) N[y|v, V ](307)

Maths for Intelligent Systems, Marc Toussaint
=N



x s S
,
y 0 0

0
0



N



101

x A⊤Q-1 a A⊤Q-1 A
,
y
Q-1 a
−Q-1 Ax s + A⊤Q-1 a S + A⊤Q-1 A,∝N
y v + Q-1 a−Q-1 A−A⊤Q-1
V + Q-1−A⊤Q-1
Q-1



N



x 0 0
,
y v 0
(308)


(309)

Entropy
H(N(a, A)) =1
log |2πeA|
2

(310)

Kullback-Leibler divergencep = N(x|a, A) ,q = N(x|b, B) ,n = dim(x) ,

 X
p(x)
D p q =
p(x) log
q(x)
x
(311)


|B|
+ tr(B -1 A) + (b − a)⊤B -1 (b − a) − n
2 D p q = log|A|
4 Dsym p q = tr(B -1 A) + tr(A-1 B) + (b − a)⊤(A-1 + B -1 )(b − a) − 2n

(312)
(313)

λ-divergence
2 Dλ p q = λ D p λp + (1−λ)q + (1−λ) D p (1−λ)p + λq

(314)

For λ = .5: Jensen-Shannon divergence.
Log-likelihoodsi1hlog|2πA| + (x-a)⊤ A-1 (x-a)
2i1hlog N[x|a, A] = − log|2πA-1 | + a⊤A-1 a + x⊤Ax − 2x⊤a2XN(x|b, B) log N(x|a, A) = −D N(b, B) N(a, A) − H(N(b, B))
log N(x|a, A) = −

(315)
(316)
(317)

x

Mixture of Gaussians
Collapsing a MoG into a single GaussianXargmin D
pi N(ai , Ai ) N(b, B)
b,B

=

(318)

i



b=

X
i

pi ai , B =Xi
pi (Ai + ai a⊤i − b b⊤)

(319)

0
V



102

B

Maths for Intelligent Systems, Marc ToussaintFurther Differential Geometry

Emphasize strong relation between a Riemannian metric (and respective geodesic)
and cost (in an optimization formulation). Pullbacks and costs. Only super brief,connections.
 Manifolds

Local tangent spaces, connection. example of kinematics
 Lie groups

exp and log
 Information Geometry[Integrate notes on information geometry]

References

J. M. Hernández-Lobato, M. W. Hoffman, and Z. Ghahramani. Predictive entropy search for efficient global optimization of black-box functions.
In
Z. Ghahramani, M. Welling, C. Cortes, N. Lawrence, and K. Weinberger,
editors, Advances in Neural Information Processing Systems 27, pages 918–
926. Curran Associates, Inc., 2014.
URL http://papers.nips.cc/paper/
5324-predictive-entropy-search-for-efficient-global-optimization-of-blackpdf.
D. Jones, M. Schonlau, and W. Welch. Efficient global optimization of expensive blackbox functions. Journal of Global Optimization, 13:455–492, 1998.
M. Toussaint. A novel augmented lagrangian approach for inequalities and convergent
any-time non-central updates. e-Print arXiv:1412.4329, 2014.

Index
L1 -norm, 6
L2 -norm, 6
Lp -norm, 6
p-norm, 6(dual problem), 62(primal problem), 62
basis, 24basis functions, 8Bayes’ Theorem, 83
Bayesian parameter estimate:, 90
Bayesian prediction: , 90
belief MDP, 71
BFGS, 55
Boltzmann distribution, 88
bound constraints, 63
characteristic polynomial, 36
column space, 32
computation graph, 10
concave, 13
conditional, 82
conjugate, 56, 83
continuous, 7
contra-variant, 40
convex, 13, 66
coordinates, 24
covariance, 85
covariance matrix, 35
covariant, 40
cumulative probability distribution, 84
damping, 52
determinant, 32
diagonalization, 36
differentiable, 7
direct sum, 6
directional derivative, 12domain, 81, 82
dual, 28
dual space, 24
dual variable, 51
dual vector, 24

eigendecomposition, 36
eigenvector, 36energy is additive, 88
Entropy is expected error, 89
error is additive, 88
Euclidean, 29
Euclidean length, 6
Expected Improvement (EI), 72
factor graph, 89
Fisher metric, 92
function network, 10
gradient, 12
Hadamard product, 7Hessian, 13
independent, 82
input null space, 32
interior point method, 60
Jacobian, 11
joint, 82
Kronecker, 6
Kronecker delta, 25
Kullback-Leibler divergence, 89
Lagrange multiplier, 63
Least Squares, 52
Levenberg-Marquardt, 52
Limited memory BFGS (L-BFGS), 55
linear, 23
Linear Program (LP):, 67
linearly independent, 24
LP in standard form:, 67
Manhattan distance, 6
marginal, 82
matrix, 26
matrix convention, 26
matrix multiplication, 27

104

Maths for Intelligent Systems, Marc ToussaintMaximum a posteriori (MAP) parameter
estimate:, 89
Maximum likelihood (ML) parameter estimate:, 89
mean, 85
metric tensor, 29
minimum for a matrix expression, 15
mixture, 86
mode, 83
modified (=approximate) KKT conditions,
65
monimials, 8
multi-linear, 23
non-parameteric, 8
orthonormal, 28
outer product, 25
output null space, 32
parameteric functions, 8particle distribution, 86
perfect line search, 56
point particle, 86
polynomial, 8polytope, 68
precision, 85
Predictive Entropy Search, 72, 73
primal-dual space, 64
probabilities are multiplicative, 88
probability density function (pdf), 84
probability distribution, 82, 84
probability measure, 81
Probability of Improvement (MPI), 72
product rule, 83
projection, 33
Quadratic Program (QP):, 67
Quadratically Constrained QP (QCQP),
67
quasiconvex, 66
random variable, 82
rank, 32
relaxations, 67
root, 49row space, 32
saddle point, 13
saddle point of the Lagrangian, 62
sample space, 81
scalar product, 28
Slater condition, 64
smooth, 7
steepest descent, 41
strong duality, 63, 64
sum-of-squares, 52
tensor product, 25
the approximate Hessian is the pullback
of a Euclidean cost feature metric, 54
total derivative, 10transpose, 30
Trust region method:, 51
unimodal, 66
Upper Confidence Bound (UCB), 72
vector, 22
vector space, 22volume, 35
Wolfe condition, 48