HT2026 Continuous Maths notes

Remaining TODOs: 27

1. Derivatives and Taylor’s Theorem

1.1. Continuity and differentiation

Definition 1.1.1: Let $𝑓 : 𝐷 \to ℝ, 𝐷 \subseteq ℝ$ .

$𝑓$ is continuous at $𝑥$ if ${lim}_{ℎ \to 0} 𝑓 (𝑥 + ℎ) = 𝑓 (𝑥)$ .

$𝑓$ is continuous if it is continuous at every $𝑥 \in 𝐷$ .

Definition 1.1.2: Let $𝑓 : 𝐷 \to ℝ, 𝐷 \subseteq ℝ^{𝑛}$ .

$𝑓$ is continuous at $𝒙$ if ${lim}_{𝒉 \to 𝟎} 𝑓 (𝒙 + 𝒉) = 𝑓 (𝒙)$ .

$𝑓$ is continuous if it is continuous at every $𝒙 \in 𝐷$ .

Theorem 1.1.3: If $𝑔$ and $𝑔$ are continuous, so are:

$𝑓 + 𝑔$
$𝑐 𝑓 \forall 𝑐 \in ℝ$
$𝑓^{𝑛} \forall 𝑛 \in ℕ$
$𝑓 𝑔$
$𝑓 \circ 𝑔$
$max (𝑓, 𝑔)$
$| 𝑓 |$
$exp (𝑓)$
$𝑓^{𝛼} \forall 𝛼 \in ℝ$ , when $𝑓$ strictly positive
$log 𝑓$ , where $𝑓$ strictly positive
$\frac{𝑓}{𝑔}$ , where $𝑔$ nonzero

Definition 1.1.4: Let $𝑓 : 𝐷 \to ℝ, 𝐷 \subseteq ℝ$ .

$𝑓$ is differentiable at $𝑥$ if ${lim}_{ℎ \to 0} \frac{𝑓 (𝑥 + ℎ) - 𝑓 (𝑥)}{ℎ} = 𝑑$ exists. Then $𝑑$ is the derivative of $𝑓$ at $𝑥$ .

$𝑓$ is differentiable if it is differentiable at every $𝑥 \in 𝐷$ .

Remark: If a function is not continuous, it is not differentiable.

Definition 1.1.5: A secant is a line that connects two points on a surface.

A tangent can be thought of as a limit of secants as the distance between the two points tend to 0.

Remark: If

𝑓

and

𝑔

are differentiable, so are all the combinations that preserve continuity mentioned earlier, with the exception of

max (𝑓, 𝑔)

and

| 𝑓 |

Remark: The derivative is a linear operator.

Theorem 1.1.6 (Quotient rule):

\frac{d}{d 𝑥} \frac{𝑓}{𝑔} = \frac{𝑔 \frac{d 𝑓}{d 𝑥} - 𝑓 \frac{d 𝑔}{d 𝑥}}{𝑔^{2}} .

Theorem 1.1.7 (Derived quotient rule):

\frac{d}{d 𝑥} (\frac{𝑓}{𝑔^{𝑘}}) = \frac{𝑔 \frac{d 𝑓}{d 𝑥} - 𝑘 𝑓 \frac{d 𝑔}{d 𝑥}}{𝑔^{𝑘 + 1}} .

Proof: Exercise; see problem sheet 0

1.2. Taylor’s theorem for univariate functions

Theorem 1.2.1 (Taylor's Theorem (univariate)): For $𝑓 : ℝ \to ℝ, 𝑥_{0} \in ℝ$ ,

𝑓 (𝑥) = \underset{Taylor polynomial}{\underset{⏟}{\sum_{𝑖 = 0}^{𝑘} \frac{{(𝑥 - 𝑥_{0})}^{𝑖}}{𝑖!} \frac{d^{𝑖} 𝑓}{d 𝑥^{𝑖}} (𝑥_{0})}} + \underset{Error term}{\underset{⏟}{\frac{{(𝑥 - 𝑥_{0})}^{𝑘 + 1}}{(𝑘 + 1)!} \frac{d^{𝑘 + 1} 𝑓}{d 𝑥^{𝑘 + 1}} (𝜉)}}

for some $𝜉 \in (𝑥_{0}, 𝑥)$ .

The Taylor polynomial (of degree $𝑘$ ) is also written ${\hat{𝑓}}_{𝑘} (𝑥)$ , and the error term, or Lagrange remainder term is denoted by $𝑒_{𝑘 + 1} (𝑥, 𝑥_{0})$ .

So, we can write $𝑓 (𝑥) = {\hat{𝑓}}_{𝑘 (𝑥)} + 𝑒_{𝑘 + 1} (𝑥, 𝑥_{0})$ for any $𝑘 \in ℕ$ so long as the first $𝑘 + 1$ derivatives of $𝑓$ all exist and are continuous on $(𝑥_{0}, 𝑥)$ .

For small $𝑒_{𝑘 + 1} (𝑥, 𝑥_{0})$ , $𝑓 (𝑥) \approx {\hat{𝑓}}_{𝑘 (𝑥)}$ .

Remark:

{\hat{𝑓}}_{𝑘} (𝑥)

is the unique

𝑘

-order polynomial that agrees with

𝑓

as to its first

𝑘

derivatives at

𝑥_{0}

Remark:Taylor’s Theorem is not the same as a Taylor series; the Taylor series is an infinite sum and does not always converge to the correct answer. (If

𝑒_{𝑘 + 1} (𝑥, 𝑥_{0})

tends to 0 for large

𝑘

and small enough

(𝑥_{0} - 𝑥)

, the series does converge and that function is analytic).

1.3. Partial and vector derivatives

Definition 1.3.1: Let $𝑓 : ℝ^{2} ⟶ ℝ$ .

$\frac{𝜕 𝑓}{𝜕 𝑥}$ is the partial derivative of $𝑓 (𝑥, 𝑦)$ with respect to $𝑥$ , as if $𝑦$ were constant.

$\frac{𝜕^{2} 𝑓}{𝜕 𝑦 𝜕 𝑥}$ means differentiate first with respect to $𝑥$ , then $𝑦$ . In general this is not equal to $\frac{𝜕^{2} 𝑓}{𝜕 𝑥 𝜕 𝑦}$ , but often it is.

Definition 1.3.2: Let $𝑓 : ℝ^{𝑛} ⟶ ℝ$ . Then the vector derivative is

\frac{d 𝑓}{d 𝒙} = (\begin{matrix} \frac{𝜕 𝑓}{𝜕 𝑥_{1}} \\ ⋮ \\ \frac{𝜕 𝑓}{𝜕 𝑥_{𝑛}} \end{matrix}) .

It turns out that this makes sense to be a derivative, because

lim_{𝒉 \to 𝟎} \frac{𝑓 (𝒙 + 𝒉) - 𝑓 (𝒙) - {(\frac{d 𝑓}{d 𝒙})}^{⊤} 𝒉}{‖ 𝒉 ‖} = 0 .

Remark: Some standard derivatives for vectors that are analogous to scalars:

$\frac{d}{d 𝒙} 𝒂^{⊤} 𝒙 = 𝒂$
$\frac{d}{d 𝒙} 𝒙^{⊤} 𝐀 𝒙 = (𝐀 + 𝐀^{⊤}) 𝒙$
$\frac{d}{d 𝒙} ‖ 𝒙 ‖ = \frac{𝒙}{‖ 𝒙 ‖}$

Remark: Similarly to scalar derivatives, the vector derivative is linear and has a product and quotient rule.

The chain rule for $𝑓 \circ 𝑔, 𝑓 : ℝ^{𝑛} ⟶ ℝ, 𝑔 : ℝ \to ℝ$ is the same as for scalars:

\frac{d}{d 𝒙} (𝑔 \circ 𝑓) = (\frac{d 𝑔}{d 𝒙} \circ 𝑓) \frac{d 𝑓}{d 𝒙} .

For $𝒇 : ℝ^{𝑚} ⟶ ℝ^{𝑛}, 𝑔 : ℝ^{𝑛} \to ℝ$ ,

\frac{d}{d 𝒙} (𝑔 \circ 𝒇) = 𝐉 {(𝒇)}^{⊤} (\frac{d 𝑔}{d 𝒙} \circ 𝒇)

where $𝐉$ is the Jacobian - see Definition 1.4.2.

Remark: The type of the vector derivative of $𝑓 : ℝ^{𝑛} \to ℝ$ is

\frac{d 𝑓}{d 𝒙} : ℝ^{𝑛} \to ℝ^{𝑛} .

Also

\frac{d}{d 𝒙} : (ℝ^{𝑛} \to ℝ) \to (ℝ^{𝑛} \to ℝ^{𝑛}) .

Definition 1.3.3: Let $𝑓 : ℝ^{𝑛} \to ℝ$ . Then the Hessian of $𝑓$ is

𝐇 (𝑓) = (\begin{matrix} \frac{𝜕^{2} 𝑓}{𝜕 𝑥_{1}^{2}} & \dots & \frac{𝜕^{2} 𝑓}{𝜕 𝑥_{1} 𝜕 𝑥_{𝑛}} \\ ⋮ & ⋱ & ⋮ \\ \frac{𝜕^{2} 𝑓}{𝜕 𝑥_{𝑛} 𝜕 𝑥_{1}} & \dots & \frac{𝜕^{2} 𝑓}{𝜕 𝑥_{𝑛}^{2}} \end{matrix}) .

Remark: The Hessian is analogous to the second derivative of a scalar function.

Theorem 1.3.4 (Clairaut's theorem): If everything is continuous,

𝐇 (𝑓) = 𝐇 {(𝑓)}^{⊤}

1.4. Vector fields and the Jacobian

Definition 1.4.1: A vector field or vector-valued function is a function

𝒇 : 𝐷 \to ℝ^{𝑚}

for

𝐷 \subseteq ℝ^{𝑛}

Remark: A vector field $𝑓$ can be decomposed into

𝒇 (𝒙) = (\begin{matrix} 𝑓_{1} (𝒙) \\ ⋮ \\ 𝑓_{𝑚} (𝒙) \end{matrix}) .

Definition 1.4.2: The derivative of a vector field $𝒇 : ℝ^{𝑛} \to ℝ^{𝑚}$ is the Jacobian, $𝐉 (𝒇) : ℝ^{𝑛} \to ℝ^{𝑚 \times 𝑛}$ , where

𝐉 (𝒇) = (\begin{matrix} \frac{𝜕 𝑓_{1}}{𝜕 𝑥_{1}} & \dots & \frac{𝜕 𝑓_{1}}{𝜕 𝑥_{𝑛}} \\ ⋮ & ⋱ & ⋮ \\ \frac{𝜕 𝑓_{𝑚}}{𝜕 𝑥_{1}} & \dots & \frac{𝜕 𝑓_{𝑚}}{𝜕 𝑥_{𝑛}} \end{matrix}) .

Remark: For $𝑓 : ℝ^{𝑛} \to ℝ$ ,

𝐉 (𝒇) = {(\frac{d 𝑓}{d 𝒙})}^{⊤} .

Remark: The Jacobian is a linear operator, i.e.

𝐉 (𝐀 𝒇) = 𝐀 𝐉 (𝒇)

for constant

𝐀

Definition 1.4.3: The Jacobian has a dot-product rule:

𝐉 (𝒇^{⊤} 𝒈) = 𝒈^{⊤} 𝐉 (𝒇) + 𝒇^{⊤} 𝐉 (𝒈) .

Definition 1.4.4: The Jacobian has a scalar-vector product rule:

𝐉 (𝑓 𝒈) = 𝒈 \frac{d 𝑓^{⊤}}{d 𝒙} + 𝑓 𝐉 (𝒈) .

Definition 1.4.5: The Jacobian has a chain rule, called the partial chain rule:

𝐉 (𝒈 \circ 𝒇) = (𝐉 (𝑔) \circ 𝒇) 𝐉 (𝒇) .

Remark: Some standard Jacobians:

\begin{matrix} 𝐉 (𝒙) = 𝐼; \\ 𝐉 (𝐀 𝒙) = 𝐀 . \end{matrix}

Remark: The Hessian is a Jacobian:

𝐇 (𝑓) = 𝐉 {(\frac{d 𝑓}{d 𝒙})}^{⊤} .

Luckily because of Clairaut’s Theorem, the transpose can usually be ignored.

Note that this basically just says that the second derivative is the derivative of the derivative.

Remark: Given vectors $𝒙_{𝑖}$ ,

\sum 𝒙_{𝑖} 𝒙_{𝑖}^{⊤} = 𝐗,

where the rows of $𝐗$ are the $𝒙_{𝑖}^{⊤}$ s.

Moreover, given scalars $𝑑_{𝑖}$ ,

\sum 𝑑_{𝑖} 𝒙_{𝑖} 𝒙_{𝑖}^{⊤} = 𝐗^{⊤} 𝐃 𝐗,

with $𝐗$ as before and $𝐃$ a diagonal matrix with the $𝑑_{𝑖}$ s as the diagonal entries.

These identities can come in handy when computing the Jacobian.

1.5. Taylor’s theorem for multivariate functions

Theorem 1.5.1 (Taylor's Theorem for multivariate functions): For $𝑓 : ℝ^{𝑛} \to ℝ, 𝑘 \in ℕ$ , if $𝒙 = 𝒙_{𝟎} + 𝒉$ ,

\begin{aligned} 𝑓 (𝒙) = \\ \begin{cases} 𝑓 (𝒙_{𝟎}) \\ + [ℎ_{1} \frac{𝜕}{𝜕 𝑥_{1}} (𝑥_{0}) + \dots + ℎ_{𝑛} \frac{𝜕}{𝜕 𝑥_{𝑛}} 𝑓] (𝒙_{𝟎}) \\ + \frac{1}{2!} [{(ℎ_{1} \frac{𝜕}{𝜕 𝑥_{1}} + \dots + ℎ_{𝑛} \frac{𝜕}{𝜕 𝑥_{𝑛}})}^{2} 𝑓] (𝒙_{𝟎}) \\ + \dots \\ + \frac{1}{𝑘!} [{(ℎ_{1} \frac{𝜕}{𝜕 𝑥_{1}} + \dots + ℎ_{𝑛} \frac{𝜕}{𝜕 𝑥_{𝑛}})}^{𝑘} 𝑓] (𝒙_{𝟎}) \end{cases}} 𝑘 th order Taylor polynomial, {\hat{𝑓}}_{𝑘} \\ + \underset{Error term}{\underset{⏟}{\frac{1}{(𝑘 + 1)!} [{(ℎ_{1} \frac{𝜕}{𝜕 𝑥_{1}} + \dots + ℎ_{𝑛} \frac{𝜕}{𝜕 𝑥_{𝑛}})}^{𝑘 + 1} 𝑓] (𝒙_{𝟎} + 𝜉 𝒉)}} for some 𝜉 \in (0, 1) . \end{aligned}

Here, the brackets with partial derivatives inside are carrying out operator algebra - they are not multiplication, but follow the same distributive rules so can be though of as so.

The Taylor polynomials up to degree 2 can be written a bit more nicely:

\begin{matrix} {\hat{𝑓}}_{0} (𝒙) = 𝑓 (𝒙_{𝟎}); \\ {\hat{𝑓}}_{1} (𝒙) = 𝑓 (𝒙_{𝟎}) + 𝒉^{⊤} \frac{d 𝑓}{d 𝒙} (𝒙_{𝟎}); \\ {\hat{𝑓}}_{2} (𝒙) = 𝑓 (𝒙_{𝟎}) + 𝒉^{⊤} \frac{d 𝑓}{d 𝒙} (𝒙_{𝟎}) + \frac{1}{2} 𝒉^{⊤} 𝐇 (𝑓) (𝒙_{𝟎}) 𝒉 . \end{matrix}

For higher powers, we need tensors which are beyond the scope of this course.

The error terms (up to $𝑘 = 2$ ) can be written in a similar fashion but evaluated at $𝒙_{𝟎} + 𝜉 𝒉$ for some $𝜉 \in (0, 1)$ .

Theorem 1.5.2: ${err}_{𝑘 + 1} (𝒙, 𝒙_{𝟎})$ is small in the sense that

lim_{‖ 𝒙 - 𝒙_{𝟎} ‖ \to 0} \frac{{err}_{𝑘 + 1} (𝒙, 𝒙_{𝟎})}{{‖ 𝒙 - 𝒙_{𝟎} ‖}^{𝑘}} = 0 .

(Proof not needed.)

2. Optimisation

Definition 2.1: Let $𝑓 : ℝ^{𝑛} \to ℝ$ .

The canonical optimisation problem is to find

min_{𝑥 \in 𝐹} 𝑓 (𝑥) or {arg min}_{𝑥 \in 𝐹} 𝑓 (𝑥)

Note that $arg min$ is a set.

$𝐹 \subset ℝ^{𝑛}$ is the feasible region.

If $𝐹 = ℝ^{𝑛}$ , the optimisation is unconstrained.

If $𝐹 \subset ℝ^{𝑛}$ , the optimisation is constrained. The standard form for constraints is:

𝐹 = {𝒙 \in ℝ^{𝑛} | \underset{equality constraints}{\underset{⏟}{𝑔_{1} (𝒙) = 0, \dots, 𝑔_{𝑙} (𝒙) = 0}}, \underset{inequality constraints}{\underset{⏟}{ℎ_{1} (𝒙) \geq 0, \dots, ℎ_{𝑚} (𝒙) \geq 0}}} .

Definition 2.2: If

𝐹

is empty, there is no solution and the optimisation problem is inconsistent.

Remark: The problem can also have no solution if, in the feasible region, the function tends to

- \infty

or there is a vertical asymptote.

2.1. Optimisation in 1 dimension

Solutions to $arg {min}_{𝑥 \in ℝ} 𝑓 (𝑥)$ satisfy the 1st-order condition $\frac{d 𝑓}{d 𝑥} = 0$ . But so do non-solutions!

Stationary points can be classified by Taylor’s theorem.

If $\frac{d 𝑓}{d 𝑥} (𝑥_{0}) = 0$ , then $𝑓 (𝑥) = 𝑓 (𝑥_{0}) + \frac{𝑥 - 𝑥_{0}}{2} \frac{d^{2} 𝑓}{d 𝑥^{2}} (𝜉)$ for some $𝜉 \in (𝑥_{0}, 𝑥)$ . So if the second derivative is positive, then for some region around $𝜉$ , $𝑓 (𝑥) > 𝑓 (𝑥_{0})$ around $𝑥_{0}$ , so $𝑥$ is a minimum. Similarly, if the second derivative is negative, then $𝑓 (𝑥) < 𝑓 (𝑥_{0})$ around $𝑥_{0}$ so $𝑥$ is a maximum.

If the second derivative is 0, this doesn’t necessarily mean that it is a point of inflection; we need to look for the first nonzero derivative. If the first nonzero derivative is odd, then it is a stationary point of inflection - since the function changes concavity around the stationary point - but if the first nonzero derivative is an even derivative, use the same rules as for the second derivative.

However, this doesn’t always work (e.g. $𝑓 (𝑥) = 𝑒^{- \frac{1}{𝑥^{2}}}, 𝑓 (0) = 0$ has all of its derivatives at 0 equal to 0).

For constrained optimisations, do the same but also consider the endpoints.

2.2. Positive/negative (semi)definiteness

Definition 2.2.1: Given a symmetric matrix $𝐀$ , $𝐀$ is

positive definite if $𝒙^{⊤} 𝐀 𝒙 > 0$ for all $𝒙 \neq 0$
positive semidefinite if $𝒙^{⊤} 𝐀 𝒙 \geq 0$ for all $𝒙 \neq 0$
negative definite if $𝒙^{⊤} 𝐀 𝒙 < 0$ for all $𝒙 \neq 0$
negative semidefinite if $𝒙^{⊤} 𝐀 𝒙 \leq 0$ for all $𝒙 \neq 0$
indefinite otherwise

Theorem 2.2.2 (Spectral theorem): If $𝐀$ is an $𝑛 \times 𝑛$ symmetric matrix, then there exists a basis (for $ℝ^{𝑛}$ ) of orthogonal eigenvectors of $𝐀$ .

Proof: See LA.

Proposition 2.2.3:

These are equivalent to certain conditions on the eigenvalues: positive definite iff all eigenvalues strictly positive, positive semidefinite iff all eigenvalues positive, etc.

Proof (for positive definiteness):

Let $𝐀$ be a symmetric matrix. By the spectral theorem, there exists an orthogonal basis of $ℝ^{𝑛}$ $⟨ 𝒙_{1}, \dots, 𝒙_{𝑛} ⟩$ ; w.l.o.g. let the $𝒙_{𝑖}$ s be unit vectors. Let $𝜆_{1}, \dots, 𝜆_{𝑛}$ be the corresponding eigenvalues.

First suppose that $𝜆_{1}, \dots, 𝜆_{𝑛} > 0$ . By the spectral theorem, $⟨ 𝒙_{1}, \dots, 𝒙_{𝑛} ⟩$ forms an orthogonal basis for $ℝ^{𝑛}$ , so for any $𝒗 \in ℝ^{𝑛}$ , $\exists 𝑐_{1}, \dots, 𝑐_{𝑛} \in ℝ$ s.t. $𝒗 = 𝑐_{1} 𝒙_{1} + \dots + 𝑐_{𝑛} 𝒙_{𝑛}$ .

Then for any $𝒗 \neq 𝟎$ ,

\begin{aligned} 𝒗^{⊤} 𝐀 𝒗 & = 𝑐_{1}^{2} 𝒙_{1}^{⊤} 𝐀 𝒙_{1} + \dots + 𝑐_{𝑛}^{2} 𝒙_{𝑛}^{⊤} 𝐀 𝒙_{𝑛} by distributivity of mat-vec mult. \\ = 𝑐_{1}^{2} 𝜆_{1} 𝒙_{1}^{⊤} 𝒙_{1} + \dots + 𝑐_{𝑛}^{2} 𝜆_{𝑛} 𝒙_{𝑛}^{⊤} 𝒙_{𝑛} \\ = 𝑐_{1}^{2} 𝜆_{1} + \dots + 𝑐_{𝑛}^{2} 𝜆_{𝑛} because the 𝒙_{𝑖} s unit vectors. \end{aligned}

Since ${(*)}^{2} \geq 0$ and the $𝜆_{𝑖}$ s are strictly positive, and there exists at least one non-zero $𝑐_{𝑖}$ , the sum of the products of positive numbers must be positive, so $𝒗^{⊤} 𝐀 𝒗$ is positive definite.

Conversely, suppose that $𝐀$ is positive definite.

Then, for all $1 \leq 𝑖 \leq 𝑛$ ,

\begin{aligned} 𝒙_{𝑖}^{⊤} 𝐀 𝒙_{𝑖} > 0 because 𝒙_{𝑖} \neq 𝟎 and 𝐀 positive definite \\ ⟹ & 𝜆_{𝑖} 𝒙_{𝑖}^{⊤} 𝒙_{𝑖} > 0 \\ ⟹ & 𝜆_{𝑖} > 0 because 𝒙_{𝑖} unit vector . \end{aligned}

Remark:

Given a symmetric $2 \times 2$ matrix $𝐴 ≔ (\begin{matrix} 𝑎 & 𝑏 \\ 𝑏 & 𝑐 \end{matrix})$ with eigenvalues $𝜆_{1}, 𝜆_{2}$ ,

𝜆_{1} 𝜆_{2} = 𝑎 𝑐 - 𝑏^{2} .

Hence

$𝑎 𝑐 - 𝑏^{2} < 0 ⟹ 𝐀$ indefinite;
$𝑎 𝑐 - 𝑏^{2} > 0 \land 𝑎 > 0 ⟹ 𝐀$ positive definite;
$𝑎 𝑐 - 𝑏^{2} = 0 \land 𝑎 + 𝑐 \geq 0 ⟹ 𝐀$ positive semidefinite
$𝑎 𝑐 - 𝑏^{2} > 0 \land 𝑎 < 0 ⟹ 𝐀$ negative definite;
$𝑎 𝑐 - 𝑏^{2} = 0 \land 𝑎 + 𝑐 \leq 0 ⟹ 𝐀$ negative semidefinite.

Definition 2.2.4: The pivots of a matrix are the entries on the diagonal when in echelon form, obtained without row multiplication or row swaps.

Remark: A symmetric matrix $𝐀$ is

positive definite if all its pivots $> 0$
positive semidefinite if all its pivots are $\geq 0$
negative definite if all its pivots $< 0$
negative semidefinite if all its pivots $\leq 0$
indefinite otherwise

but the definiteness conditions hold only if the pivots can be found without row swaps

Remark: Definiteness works a bit like signs for scalars. Pos def equivalent to positive, pos semidef quivalent to nonnegative, etc. Incl. (pos def)^⁻¹ pos def (and this always exists, in this case).

Remark: If $𝐀$ is positive (semi)definite then so are

$𝑐 𝐀$ for any $𝑐 > 0$
$𝐀^{- 1}$ , which is guaranteed to exist if $𝐀$ is positive definite
$𝐀 𝐁 𝐀$ for any positive (semi)definite $𝐵$
any upper-left submatrix of $𝐀$ TODO intuition
$𝐂^{⊤} 𝐀 𝐂$ , if $𝐂$ is of full rank, regardless of if $𝐂$ is square or not - but if $𝑐$ not of full rank, we can still guarantee positive semidefiniteness TODO intuition
$𝐀 - 𝒗 𝒗^{⊤}$ , if $𝒗^{⊤} 𝐀^{- 1} 𝒗 < 1$

Remark: For any symmetric matrix

𝐂

𝐂 - 𝜆 𝐈

is positive semidefinite if

𝜆 \leq

the smallest eigenvalue of

𝐂

. It is positive definite if the inequality is strict.

Remark: The zero matrix is both positive semidefinite and negative semidefinite.

TODO: lookup: diagonally dominant

2.3. Unconstrained optimisation over $ℝ^{𝑛}$

Solutions to $\underset{𝑥 \in ℝ^{𝑛}}{arg min}$ , if they exist, satisfy the 1st-order condition $\frac{d 𝑓}{d 𝒙} = 𝟎$ . This is in general not easy to solve as it is a system of $𝑛$ (not necessarily linear) equations.

In higher dimensions, stationary points can be local minima, local maxima or saddle points.

Definition 2.3.1: A saddle point is a stationary point at which moving in some directions gives a more positive result, and in another direction decreases.

We can classify these higher-dimensional stationary points using the Hessian, similarly to how we used the second derivative for scalar functions.

Theorem 2.3.2:If

𝐇 (𝑓) (𝒙)

is positive definite, then

𝐇 (𝑓) (𝒙^{*})

is positive definite for

𝒙^{*}

close to

𝒙

Hence if the Hessian at $𝒙$ is positive definite, the stationary point at $𝒙$ is a local minimum. If the Hessian is negative definite, the stationary point is a local maximum. If the Hessian is indefinite, the stationary point is a saddle point.

If the Hessian is semidefinite, you cannot determine anything other than that the point is not a local minimum/maximum. For better analysis we would need to use the 3rd-order Taylor approximation.

Fact: the outer product of a vector with itself is always positive semidefinite. Proof: $𝒗^{⊤} (𝒃 𝒃^{⊤}) 𝒗 = {(𝒗^{⊤} 𝒃)}^{2} \geq 0$ .

2.4. Convexity

Definition 2.4.1: A set $𝐷 \subseteq ℝ^{𝑛}$ is convex if

\begin{matrix} \forall 𝒙_{𝟏}, 𝒙_{𝟐} \in 𝐷, 𝛼 \in (0, 1), \\ (1 - 𝛼) 𝒙_{𝟏} + 𝛼 𝒙_{𝟐} \in 𝐷 . \end{matrix}

That is, the line connecting any $𝒙_{𝟏}, 𝒙_{𝟐} \in 𝐷$ lies wholly in $𝐷$ .

Definition 2.4.2: A function $𝑓 : 𝐷 ⟶ ℝ$ for convex $𝐷 \subseteq ℝ^{𝑛}$ is convex if

\begin{matrix} \forall 𝒙_{𝟏}, 𝒙_{𝟐} \in 𝐷, 𝛼 \in (0, 1), \\ 𝑓 ((1 - 𝛼) 𝒙_{𝟏} + 𝛼 𝒙_{𝟐}) \leq (1 - 𝛼) 𝑓 (𝒙_{𝟏}) + 𝛼 𝑓 (𝒙_{𝟐}) . \end{matrix}

That is, any point on $𝑓$ between $𝒙_{𝟏}$ and $𝒙_{𝟐}$ lies below the secant between $𝒙_{𝟏}$ and $𝒙_{𝟐}$ .

Equivalently, every secant on $𝑓$ lies above the surface of $𝑓$ .

If the inequality is strict, $𝑓$ is strictly convex if the inequality is the other way around, $𝑓$ is concave or strictly concave in a similar manner.

Theorem 2.4.3: If

𝑓

is convex, every stationary point is a global minimum.

Theorem 2.4.4: If

𝑓

is strictly convex, there is at most one stationary point, and that is the global minimum.

Remark: If carrying out unconstrained optimisation on a convex function, just solve the first-order condition (once you have established that there is a minimum).

Remark: In one dimension,

𝑓

is convex if

\frac{d^{2} 𝑓}{d 𝑓^{2}} \geq 0

everywhere; if the inequality is strict then

𝑓

is strictly convex. The implication does not hold in reverse.

Remark: Convexity does not imply the existence of a minimum!

Remark: For higher dimensions,

𝑓

is is convex iff

𝐇 (𝑓)

is positive semidefinite everywhere; if (but not only if)

𝐇 (𝑓)

is positive definite everywhere, then

𝑓

is strictly convex.

Remark:

Some common convex and concave functions:

Linear functions are both convex and concave.
${| 𝑥 |}^{𝑝}$ for $𝑝 > 1$ is strictly convex
$𝑒^{𝑥}$ is strictly convex
$log 𝑥$ is strictly concave
for convex $𝑓$ , $- 𝑓$ is concave
$𝒙 \mapsto 𝒙^{⊤} 𝐀 𝒙$ is convex if $𝐀$ is positive semidefinite, strictly convex if $𝐀$ is positive definite
for convex $𝑓, 𝑔$ , $𝑓 + 𝑔$ is convex
for convex $𝑓$ , $𝑓 (𝐀 𝒙 + 𝒃)$ is convex (but strict convexity is only preserved if $𝐀$ has full rank)
for convex $𝑓, 𝑔$ , $max (𝑓, 𝑔)$ is convex
for convex $𝑓$ , $exp (𝑓)$ is convex

Theorem 2.4.5: For $𝑓 : ℝ ⟶ ℝ, 𝑔 : ℝ \to ℝ$ ,

\begin{matrix} 𝑓 convex \land {\begin{matrix} 𝑓 increasing \land 𝑔 convex \\ 𝑓 decreasing \land 𝑔 concave \end{matrix}} ⟹ 𝑓 \circ 𝑔 convex; \\ 𝑓 concave \land {\begin{matrix} 𝑓 increasing \land 𝑔 concave \\ 𝑓 decreasing \land 𝑔 convex \end{matrix}} ⟹ 𝑓 \circ 𝑔 concave . \end{matrix}

Proof: $𝑓 \circ 𝑔$ convex iff ${(𝑓 \circ 𝑔)}^{″} \geq 0$ everywhere.

\begin{matrix} \frac{d}{d 𝑥} (𝑓 \circ 𝑔) = (\frac{d 𝑓}{d 𝑥} (𝑔 (𝑥))) \cdot \frac{d 𝑔}{d 𝑥} . \\ \frac{d^{2}}{d 𝑥^{2}} (𝑓 \circ 𝑔) = (\frac{d^{2} 𝑓}{d 𝑥^{2}} (𝑔 (𝑥))) \cdot {(\frac{d 𝑔}{d 𝑥})}^{2} + \frac{d 𝑓}{d 𝑥} (𝑔 (𝑥)) \frac{d^{2} 𝑔}{d 𝑔^{2}} . \end{matrix}

${(\frac{d 𝑔}{d 𝑥})}^{2}$ is always positive, so we require $\frac{d^{2} 𝑓}{d 𝑥^{2}}$ to be positive everywhere, i.e. $𝑓$ is convex. We then require $\frac{d 𝑓}{d 𝑥}$ and $\frac{d^{2} 𝑓}{d 𝑔^{2}}$ to have the same signs everywhere, so either $𝑓$ is increasing and $𝑔$ is convex, or $𝑓$ is decreasing and $𝑔$ is concave.

The proof for concavity is similar.

Remark: These conditions are sufficient but not necessary.

Remark: Theorem 2.4.5 can be used also for

𝑓 : ℝ \to ℝ, 𝑔 : ℝ^{𝑛} \to ℝ

Theorem 2.4.6:

Let $𝑓 : ℝ^{𝑛} \to ℝ$ , $𝒈 : ℝ^{𝑚} \to ℝ^{𝑛}$ such that

𝒈 ≔ (\begin{matrix} 𝑔_{1} \\ ⋮ \\ 𝑔_{𝑛} \end{matrix})

with $𝑔_{1}, \dots, 𝑔_{𝑛} : ℝ^{𝑚} \to ℝ$ .

Then if $𝑓$ is (strictly) convex, and for each $𝑖 \in {1, \dots, 𝑛}$ , either

{\begin{cases} 𝑓 (strictly) increasing in its 𝑖 th argument \land 𝑔_{𝑖} (strictly) convex, or \\ 𝑓 (strictly) decreasing in its 𝑖 th argument \land 𝑔_{𝑖} (strictly) concave, \end{cases}

then $𝑓 \circ 𝒈$ is (strictly) convex.

Theorem 2.4.7 (Jenson's inequality): For a random variable $𝑋$ and convex function $𝑓$ ,

𝑓 (𝔼 [𝑋]) \leq 𝔼 [𝑓 (𝑥)] .

(This is not on the syllabus but is useful)

2.5. Optimisation tricks

Theorem 2.5.1:If $𝑔$ is strictly increasing, then

arg min 𝑓 (𝑥) = arg min 𝑔 (𝑓 (𝑥)) .

Theorem 2.5.2: If $𝑔$ is injective/1-to-1 then

arg min 𝑓 (𝑥) = 𝑔 (arg min 𝑓 (𝑔 (𝑥))) .

2.6. Optimisation over $ℝ^{𝑛}$ with equality constraints

Remark: Given a minimisation problem with an equality constraint, it can sometimes be reduced to an unconstrained problem:

Given $𝑔 (𝑥, 𝑦) = 𝑥 + 𝑦 + 1$ ,
$min_{𝒙 \in ℝ^{2} | 𝑔 (𝒙) = 0} 𝑓 (𝒙) = min_{𝑥 \in ℝ} 𝑓 (𝑥, - 1 - 𝑥)$
Given $𝑔 (𝑥, 𝑦) = 𝑥^{2} + 𝑦^{2} - 1$ ,
$min_{𝒙 \in ℝ^{2} | 𝑔 (𝒙) = 0} 𝑓 (𝑥_{1}, 𝑥_{2}) = min_{𝜃 \in (0, 2 𝜋]} 𝑓 (sin 𝜃, cos 𝜃)$

But usually, this isn’t the case.

Theorem 2.6.1 (Lagrange's theorem / method of Lagrange multipliers): For a minimisation problem with one equality constraint, minimise $𝑓 (𝒙)$ subject to $𝑔 (𝒙) = 0$ , the stationary points must satisfy

\frac{d 𝑓}{d 𝒙} = 𝜆 \frac{d 𝑔}{d 𝒙},

i.e. a first order condition is $\frac{d}{d 𝑥} (𝑓 (𝒙) - 𝜆 𝑔 (𝒙)) = 𝟎$ .

Proof:

TODO

Definition 2.6.2: The Lagrangian of such a minimisation problem is

Λ (𝜆, 𝒙) = 𝑓 (𝒙) - 𝜆 𝑔 (𝒙) .

The stationary points of $𝑓 (𝒙)$ where $𝑔 (𝒙) = 𝟎$ are the stationary points of $Λ$ .

Remark: It is not simple to determine if the stationary points of the Lagrangian are a minimum, maximum or saddle point of the problem. It is not sufficient to look at the definiteness of the Hessian of

Λ

, or the convexity of

𝑓

TODO: look up bordered Hessian

Theorem 2.6.3: For $min 𝑓 (𝒙) s.t. 𝑔_{1} (𝒙) = 0, \dots, 𝑔_{𝑙} (𝒙) = 0$ ,

the Lagrangian is

Λ (𝝀, 𝒙) = 𝑓 (𝒙) - 𝜆_{1} 𝑔_{1} (𝒙) - \dots - 𝜆_{𝑙} 𝑔_{𝑙} (𝒙) .

Then the stationary points of the minimisation problem are the stationary points of the Lagrangian.

2.7. Inequality contrained optimisation (over $ℝ^{𝑛}$ )

Definition 2.7.1: Given a minimisation problem $min 𝑓 (𝒙)$ subject to $ℎ (𝒙) \geq 0$ , at the solution, $ℎ$ is tight if $ℎ (𝒙) = 0$ , and slack if $ℎ (𝒙) > 0$ .

To solve such a minimisation, solve:

for tight constraints via the previous section;
for slack constraints, by solving the usual first-order condition $\frac{d 𝑓}{d 𝒙} = 𝟎$ (as usual), and ignoring any solutions outside of the feasible region.

Either way, $\frac{d 𝑓}{d 𝒙} = 𝜇 \frac{d ℎ}{d 𝒙}$ , with $𝜇 \geq 0$ : if we find a minimum of $𝑓$ on the boundary of the feasible region, where $ℎ (𝒙) = 0$ , then the gradients of $𝑓$ and $ℎ$ at that point must be parallel and in the same direction, otherwise we could move into the interior of the feasible region, where $ℎ (𝒙) \geq 0$ , while decreasing $𝑓$ (so on the boundary $𝜇 > 0$ ); on the interior of the feasible region, $𝑓$ can attain a minimum independent of the constraint, so $𝜇 = 0$ .

Note that in both cases, $𝜇 ℎ (𝒙) = 0$ .

Therefore for inequality constraints, we can find stationary points of the lagrangian $Λ (𝜇, 𝒙) = 𝑓 (𝒙) - 𝜇 ℎ (𝒙)$ that satisfy the constraint, simultaneously with $𝜇 \geq 0$ , $𝜇 ℎ (𝒙) = 0$ .

Corollary 2.7.2 (Complementary slackness): At a minimum of an minimisation problem with constraint $ℎ (𝒙) = 0$ and Lagrange multiplier $𝜇$ ,

[𝜇 \geq 0 \land ℎ (𝒙) = 0] \lor [𝜇 = 0 \land ℎ (𝒙) > 0] .

Theorem 2.7.3: When optimising

𝑓

with one inequality constraint, if

𝑓

is convex and the feasible set is convex, then either every stationary point inside the constraint is a global minimum, or one of the stationary points on the constraint is a global minimum.

TODO KKT conditions (not on syllabus)

TODO: what to do when mix of equality and inequality constraints

3. Algorithms for numerical integration

3.1. 1-dimensional integration

Here we will consider approximating the integral of $𝑓$ over a small interval $[𝑎, 𝑏]$ , with $𝑚 = \frac{𝑎 + 𝑏}{2}$ .

The simplest approximation is the 0th order Taylor approximation; this is equal to the approximation using the 1st order Taylor approximation.

Definition 3.1.1: The midpoint rule: define

𝑀_{1} [𝑓, 𝑎, 𝑏] ≔ \int_{𝑎}^{𝑏} {\hat{𝑓}}_{1} (𝑥) d 𝑥 = (𝑏 - 𝑎) 𝑓 (𝑚) .

The subscript 1 indicates that it is operating over 1 strip. The notation for this (and other methods in this section) is not standard.

Theorem 3.1.2: The error of the midpoint rule,

err (𝑀_{1}) [𝑓, 𝑎, 𝑏] = 𝑀_{1} [𝑓, 𝑎, 𝑏] - \int_{𝑎}^{𝑏} 𝑓 (𝑥) d 𝑥,

is bounded by

- \frac{{(𝑏 - 𝑎)}^{3}}{24} 𝐷_{2} \leq err (𝑀_{1}) [𝑓, 𝑎, 𝑏] \leq - \frac{{(𝑏 - 𝑎)}^{3}}{24} 𝐷_{2},

where

\begin{matrix} 𝐷_{𝑘} = min_{𝑥 \in (𝑎, 𝑏)} \frac{d^{𝑘} 𝑓}{d 𝑥^{𝑘}}, \\ 𝐷_{𝑘} = max_{𝑥 \in (𝑎, 𝑏)} \frac{d^{𝑘} 𝑓}{d 𝑥^{𝑘}} . \end{matrix}

Proof:

\begin{aligned} err (𝑀_{1}) [𝑓, 𝑎, 𝑏] & = 𝑀_{1} [𝑓, 𝑎, 𝑏] - \int_{𝑎}^{𝑏} 𝑓 (𝑥) d 𝑥 \\ = (𝑏 - 𝑎) 𝑓 (𝑚) - \int_{𝑎}^{𝑏} 𝑓 (𝑚) + (𝑥 - 𝑚) \frac{d 𝑓}{d 𝑥} (𝑚) + \frac{{(𝑥 - 𝑚)}^{2}}{2} \frac{d^{2} 𝑓}{d 𝑥^{2}} (𝜉) d 𝑥 \\ for some 𝜉 \in (𝑎, 𝑏) \\ = (𝑏 - 𝑎) 𝑓 (𝑚) - (𝑏 - 𝑎) 𝑓 (𝑚) - {[(𝑥 - 𝑚)]}_{𝑎}^{𝑏} \frac{d 𝑓}{d 𝑥} (𝑚) - {[\frac{{(𝑥 - 𝑚)}^{3}}{6}]}_{𝑎}^{𝑏} \frac{d^{2} 𝑓}{d 𝑥^{2}} (𝜉) \\ = - (\frac{{(𝑏 - \frac{𝑎 + 𝑏}{2})}^{3}}{6} - \frac{{(𝑎 - \frac{𝑎 + 𝑏}{2})}^{3}}{6}) \frac{d^{2} 𝑓}{d 𝑥^{2}} (𝜉) \\ = - \frac{{(𝑏 - 𝑎)}^{3}}{24} \frac{d^{2} 𝑓}{d 𝑥^{2}} (𝜉) . \end{aligned}

We can bound $\frac{d^{2} 𝑓}{d 𝑥^{2}} (𝜉)$ by $𝐷_{2} \leq \frac{d^{2} 𝑓}{d 𝑥^{2}} (𝜉) \leq 𝐷_{2}$ , so

- \frac{{(𝑏 - 𝑎)}^{3}}{24} 𝐷_{2} \leq err (𝑀_{1}) [𝑓, 𝑎, 𝑏] \leq - \frac{{(𝑏 - 𝑎)}^{3}}{24} 𝐷_{2} .

Definition 3.1.3: Instead of Taylor’s theorem, we can use polynomial interpolation - a polynomial that agrees with

𝑓

at some points.

Definition 3.1.4: The Trapezium rule creates a trapezium with the axis and the endpoints of the strip. It is generally worse than the midpoint rule, so will not be assessed in the course.

𝑇_{1} [𝑓, 𝑎, 𝑏] = \frac{𝑏 - 𝑎}{2} (𝑓 (𝑎) + 𝑓 (𝑏)) .

Theorem 3.1.5: The error bound for the trapezium rule is

\frac{1}{12} {(𝑏 - 𝑎)}^{3} 𝐷_{2} \leq err (𝑇_{1}) [𝑓, 𝑎, 𝑏] \leq \frac{1}{12} {(𝑏 - 𝑎)}^{3} 𝐷_{2} .

Remark: The error bound is twice as big as the midpoint rule.

Definition 3.1.6: Simpson’s rule interpolates a parabola from the endpoints and midpoints. The derivation is quite complicated but the resulting formula is very simple: it is just a weighted average. Although this is over a single strip, Simpson’s rule is usually considered to operate over 2 strips, for notational consistency (as we will see later).

𝑆_{2} [𝑓, 𝑎, 𝑏] = \frac{𝑏 - 𝑎}{6} (𝑓 (𝑎) + 4 𝑓 (𝑚) + 𝑓 (𝑏)) .

Theorem 3.1.7: The error bound for Simpson’s rule is

\frac{1}{2880} {(𝑏 - 𝑎)}^{5} 𝐷_{4} \leq err (𝑆_{2}) [𝑓, 𝑎, 𝑏] \leq \frac{1}{2880} {(𝑏 - 𝑎)}^{5} 𝐷_{4} .

Remark: This error bound is much better, if

𝑎

and

𝑏

are close. But this error bound only applies if the function has at least 4 derivatives; otherwise, the method can still be used, but we can’t bound the error in this manner.

Remark:Boole’s rule gives a similar formula for interpolating a quartic polynomial, but is a bit more complicated and doesn’t offer much improvement over Simpson’s rule.

Higher order Taylor approximations might not be a great choice for a better approximation because:

you need to evaluate high derivatives
Taylor approximations are local
functions might not have enough derivatives/be analytic

Higher-order polynomial interpolation also might not be a good idea:

it is not necessarily numerically unstable (as commonly believed), but it can be if you do it badly
the Runge phenomenon - some functions get worse and worse approximations as the order of the interpolated polynomial increases because polynomial interpolation works worse near the ends of the interval
but Chebyshev interpolations are much better! These have interpolation points closer together near to the ends of the interval

Definition 3.1.8: The Runge function is

𝑓 (𝑥) = \frac{1}{1 + 25 {(2 𝑥 - 1)}^{2}} .

Higher-order polynomial interpolation on the Runge function, using evenly spaced points, is really bad.

A better method is just to split the interval into smaller strips and use an easier method on each of them. This is guaranteed to converge to the right answer regardless of which method you use, so long as the function is continuous.

Definition 3.1.9: The composite midpoint rule is essential just an average of all the heights of the $𝑛 + 1$ endpoints of $𝑛$ strips, multiplied by the width of the region.

𝑀_{𝑛} [𝑓, 𝑎, 𝑏] = \frac{𝑏 - 𝑎}{𝑛} (𝑓 (\frac{𝑥_{0} + 𝑥_{1}}{2}) + \dots + 𝑓 (\frac{𝑥_{𝑛 - 1} + 𝑥_{𝑛}}{2})) .

Theorem 3.1.10: The error bound for the composite midpoint rule is

- \frac{{(𝑏 - 𝑎)}^{3}}{24 𝑛^{2}} 𝐷_{2} \leq err (𝑀_{𝑛}) [𝑓, 𝑎, 𝑏] \leq - \frac{{(𝑏 - 𝑎)}^{3}}{24 𝑛^{2}} 𝐷_{2} .

Remark: The main computation cost is

𝑛

evaluations of

𝑓

; the error is

𝑂 (𝑛^{- 2})

so this is pretty good for the amount of work we do.

Definition 3.1.11: The composite Simpson’s rule is again akin to a weighted average, now with coefficients $1, 4, 2, 4, \dots, 4, 1$ , with the 1s and 4s as before and the 2s coming from where the endpoints are repeated.

𝑆_{𝑛} [𝑓, 𝑎, 𝑏] = \frac{𝑏 - 𝑎}{3 𝑛} (𝑓 (𝑥_{0}) + 4 𝑓 (𝑥_{1}) + 2 𝑓 (𝑥_{2}) + 4 𝑓 (𝑥_{3}) + \dots + 4 𝑓 (𝑥_{𝑛 - 1}) + 𝑓 (𝑥_{𝑛})) .

Theorem 3.1.12: The error bound for the composite Simpson’s rule is

\frac{{(𝑏 - 𝑎)}^{5}}{180 𝑛^{4}} 𝐷_{4} \leq err (𝑆_{𝑛}) [𝑓, 𝑎, 𝑏] \leq \frac{{(𝑏 - 𝑎)}^{5}}{180 𝑛^{4}} 𝐷_{4} .

Remark: The main computational cost is again in evaluations of

𝑓

, this time

𝑛 + 1 = 𝑂 (𝑛)

. But now the error bound is

𝑂 (𝑛^{- 4})

which is even better relative to the amount of work we do.

We can make improvements on these methods:

Pick more optimal strip endpoints (Gaussian quadrature - TODO what is this?)
Extrapolate a sequence of $𝑇_{𝑛}$ or $𝑀_{𝑛}$ (Romberg integration - TODO what is this?)
Recursive subdivision (reuse evaluation, do more [via thinner strips] on difficult bits) - based on 4th derivative being large, or use the a posteriori error estimate
$| err (𝑆_{2 𝑛}) [𝑓, 𝑎, 𝑏] | \approx \frac{1}{16} | err (𝑆_{𝑛}) [𝑓, 𝑎, 𝑏] | \approx \frac{1}{15} | 𝑆_{2 𝑛} [𝑓, 𝑎, 𝑏] | - 𝑆_{𝑛} [𝑓, 𝑎, 𝑏]$
TODO proof, how, why??
Or just do more work where $𝑓$ is large, as the errors have a greater effect there!

Also note that these methods can’t help with integrals with $\pm \infty$ in either of the limits.

3.2. Integration in $𝑑$ dimensions

In two dimensions, $\int_{𝑅} 𝑓 (𝑥, 𝑦) d (𝑥, 𝑦)$ is the volume under a surface in region $𝑅 \subseteq ℝ^{2}$ .

Theorem 3.2.1 (Fubini's Theorem):

\int_{𝑅} 𝑓 (𝑥, 𝑦) d (𝑥, 𝑦) = \int_{𝑥_{0}}^{𝑥_{1}} \int_{𝑦_{0}}^{𝑦_{1}} 𝑓 (𝑥, 𝑦) d 𝑦 d 𝑥

for sufficiently nice $𝑓$ and $𝑅$ (we need to be able to express the bounds nicely). This “split” integral is called an iterated integral.

This extends to higher dimensions as expected.

As long as everything is integrable, the order doesn’t matter.

We can use the methods we already know to approximate higher-dimensional integrals, as long as we can split the integral into iterated integrals.

Remark: The midpoint rule can be extended to multiple dimensions by converting into an iterated integral via Fubini’s theorem, and then applying the 1-dimensional midpoint rule to successive inner integrals. For a $𝑑$ -dimensional integral $𝐼$ of the function $𝑓$ over the unit hypercube, the midpoint rule with $𝑛$ strips in each dimension gives

𝑀_{𝑛} [𝑓, 𝟎, 𝟏] = \frac{1}{𝑛^{𝑑}} \sum_{𝑖_{1} = 0}^{𝑛} \dots \sum_{𝑖_{𝑑} = 0}^{𝑛} 𝑓 (\frac{2 𝑖_{1} - 1}{2 𝑛}, \dots, \frac{2 𝑖_{𝑑} - 1}{2 𝑛}) .

This is just taking the average of the value of $𝑓$ at the midpoints of each hypercube in a grid over the region.

Remark: Simpson’s rule can also be extended to higher dimensions: for 2 dimensions,

𝑆_{𝑛} [𝑓, 𝟎, 𝟏] = \frac{1}{9 𝑛^{2}} \sum_{𝑖 = 0}^{𝑛 + 1} \sum_{𝑗 = 0}^{𝑛 + 1} 𝑤_{𝑖 𝑗} 𝑓 (\frac{𝑖}{𝑛 + 1}, \frac{𝑗}{𝑛 + 1})

where

𝐖 = (\begin{matrix} 1 & 4 & 2 & 4 & \dots & 4 & 1 \\ 4 & 16 & 8 & 16 & \dots & 16 & 4 \\ 2 & 8 & 4 & 8 & \dots & 8 & 2 \\ 4 & 16 & 8 & 16 & \dots & 16 & 4 \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋱ & ⋮ & ⋮ \\ 4 & 16 & 8 & 16 & \dots & 16 & 4 \\ 1 & 4 & 2 & 4 & \dots & 4 & 1 \end{matrix}) .

Remark: In $𝑑$ dimensions,

\begin{matrix} err (𝑀_{𝑛}) [𝑓, 𝑅] = 𝑂 (𝑑 𝑛^{- 2}) \\ err (𝑆_{𝑛}) [𝑓, 𝑅] = 𝑂 (𝑑 𝑛^{- 4}) . \end{matrix}

This seems quite good, but we do $𝑂 (𝑛^{𝑑})$ evaluations of the function; so relative to the work, the error becomes steadily worse.

This is called the “curse of dimensionality”… all of these methods improve extremely slowly in high dimensions.

Definition 3.2.2: Monte Carlo integration does the midpoint rule, but using random points according to some distribution.

Consider sampling random points along a scalar function:

\begin{matrix} 𝑋 \sim Unif [𝑎, 𝑏] (i.e. p.d.f 𝑝 (𝑥) = \frac{1}{𝑏 - 𝑎} for 𝑥 \in [𝑎, 𝑏], 0 otherwise) . \\ Note that 𝔼 [𝑓 (𝑋)] = \int_{- \infty}^{\infty} 𝑓 (𝑥) 𝑝 (𝑥) d 𝑥 = \frac{1}{𝑏 - 𝑎} \int_{𝑎}^{𝑏} 𝑓 (𝑥) d 𝑥, \\ and also that lim_{𝑁 \to \infty} [\frac{1}{𝑁} \sum_{𝑖 = 1}^{𝑁} 𝑓 (𝑥_{𝑖})] = 𝔼 [𝑓 (𝑋)] . \end{matrix}

Therefore, given $𝑋_{1}, \dots, 𝑋_{𝑁}$ i.i.d. according to $Unif [𝑎, 𝑏]$ , for large enough $𝑁$ ,

\int_{𝑎}^{𝑏} 𝑓 (𝑥) d 𝑥 \approx \frac{𝑏 - 𝑎}{𝑁} \sum_{𝑖 = 1}^{𝑁} 𝑓 (𝑋_{𝑖}) .

We extend this to $𝑑$ dimensions, using the notation

𝑀 𝐶_{𝑁} [𝑓, 𝑅] ≔ 𝐴 (𝑅) \frac{1}{𝑁} \sum_{𝑖 = 1}^{𝑁} 𝑓 (𝑿_{𝒊})

for large $𝑁$ to approximate $\int_{𝑅} 𝑓 (𝒙) d 𝒙$ , with $𝑿_{1}, \dots, 𝑿_{𝑁}$ i.i.d. uniformly distributed across the region $𝑅$ , and $𝐴 (𝑅)$ the area of the region $𝑅$ .

Remark: This is called Monte Carlo integration because it is a Monte Carlo algorithm; the answer is random.

Remark: This is in contrast to Las Vegas algorithms, which have random running time but nonrandom answer.

Lemma 3.2.3:

Monte Carlo integration is unbiased (on average it is correct).

Proof: Take $𝑿_{𝟏}, \dots, 𝑿_{𝑵}$ i.i.d. with pdf

𝑝 (𝒙) = {\begin{cases} \frac{1}{𝐴 (𝑅)} & for 𝒙 \in 𝑅, \\ 0 & otherwise . \end{cases}

Then

\begin{aligned} 𝔼 [𝑀 𝐶_{𝑁} [𝑓, 𝑅]] & = 𝔼 [𝐴 (𝑅) \frac{1}{𝑁} \sum_{𝑖 = 1}^{𝑁} 𝑓 (𝑿_{𝒊})] \\ = \frac{𝐴 (𝑅)}{𝑁} \sum_{𝑖 = 1}^{𝑁} 𝔼 [𝑓 (𝑿_{𝒊})] by linearity of expectation \\ = \frac{𝐴 (𝑅)}{𝑁} \sum_{𝑖 = 1}^{𝑁} \int_{𝑅} \frac{𝑓 (𝒙)}{𝐴 (𝑅)} d 𝒙 \\ = \int_{𝑅} 𝑓 (𝒙) d 𝒙 . \end{aligned}

Proposition 3.2.4:

Var (𝑀 𝐶_{𝑁} [𝑓, 𝑅]) = \frac{𝑉}{𝑁}

for a constant $𝑉$ . Moreover, $𝑉$ can be approximated by $\hat{𝑉}$ that can be expressed in terms of $𝐴 (𝑅)$ , $𝑁$ , and the $𝑓 (𝑿_{𝒊})$ s.

Proof:

\begin{aligned} Var (𝑀 𝐶_{𝑁} [𝑓, 𝑅]) & = Var (𝐴 (𝑅) \frac{1}{𝑁} \sum 𝑓 (𝑿_{𝒊})) \\ = \frac{{𝐴 (𝑅)}^{2}}{𝑁^{2}} Var (\sum 𝑓 (𝑿_{𝒊})) \\ = \frac{{𝐴 (𝑅)}^{2}}{𝑁^{2}} \sum Var (𝑓 (𝑿_{𝒊})) because 𝑋_{𝑖} s indep because 𝑿_{𝒊} s i.i.d. \\ = \frac{{𝐴 (𝑅)}^{2}}{𝑁} Var (𝑓 (𝑿_{𝒊})) because all variances equal because 𝑿_{𝒊} s i.i.d. \\ = \frac{𝑉}{𝑁} with 𝑉 ≔ {𝐴 (𝑅)}^{2} Var (𝑓 (𝑿_{𝒊})) . \end{aligned}

\begin{aligned} 𝑉 & = {𝐴 (𝑅)}^{2} Var (𝑓 (𝒙)) \\ = {𝐴 (𝑅)}^{2} (𝔼 [{𝑓 (𝑿_{𝒊})}^{2}] - {𝔼 [𝑓 (𝑿_{𝒊})]}^{2}) \\ = {𝐴 (𝑅)}^{2} (\int_{𝑅} \frac{{𝑓 (𝒙)}^{2}}{𝐴 (𝑅)} d 𝒙 - {(\int_{𝑅} \frac{𝑓 (𝒙)}{𝐴 (𝑅)} d 𝒙)}^{2}) \\ \approx \frac{{𝐴 (𝑅)}^{2}}{𝑁} \sum_{𝑖} 𝑓^{2} (𝑿_{𝒊}) - {(\frac{𝐴 (𝑅)}{𝑁} \sum_{𝑖} 𝑓 (𝑿_{𝒊}))}^{2} \\ ≕ \hat{𝑉} . \end{aligned}

Definition 3.2.5: The standard error is the square root of the estimated variance,

SE = \sqrt{\frac{\hat{𝑉}}{𝑁}} = 𝑂 (𝑁^{- \frac{1}{2}})

Remark: This all looks quite bad, but none of it is dependent on the dimension, so this is actually quite good! The goodness of the algorithm only depends on the area of the region we are integrating over, and the “wiggliness” of the function.

Lemma 3.2.6:

Monte Carlo integration is consistent:

ℙ (| err (𝑀 𝐶_{𝑁}) [𝑓, 𝑅] | \geq 𝜀) \underset{𝑁 \to \infty}{\to} 0 \forall 𝜀 .

Proof: Let $𝑌 = 𝑀 𝐶_{𝑁} [𝑓, 𝑅]$ .

Then

𝔼 [𝑌] = \int_{𝑅} 𝑓 (𝒙) d 𝒙, Var (𝑌) = \frac{𝑉}{𝑁} for some 𝑉 \in ℝ .

Then by Chebyshev’s inequality,

ℙ (| 𝑌 - 𝔼 [𝑌] | \geq 𝑐) \leq \frac{𝑉}{𝑁 𝑐^{2}} \underset{𝑁 \to \infty}{\to} 0

Theorem 3.2.7 (Central Limit Theorem): If $𝑌_{1}, . ., 𝑌_{𝑁}$ i.i.d., then

𝔼 [𝑌] = 𝜇 \land Var (𝑌) = 𝜎^{2} > 0 ⟹ \frac{1}{𝑁} \sum_{𝑖 = 1}^{𝑁} 𝑌_{𝑖} \to 𝑁 (𝜇, \frac{𝜎^{2}}{𝑁})

(This is not a precise statement.)

That is, a large sample of i.i.d. random variables tends to a normal distribution.

Lemma 3.2.8:

For large $𝑁$ , the error is distributed approximately according to $𝑁 (0, \frac{\hat{𝑉}}{𝑁})$ .

Proof:

\begin{aligned} Let 𝑌_{𝑖} & = 𝑓 (𝑿_{𝒊}) 𝐴 (𝑅) - \int_{𝑅} 𝑓 (𝒙) d 𝒙 \\ = err (𝑀 𝐶_{1}) [𝑓, 𝑅] . \end{aligned}

Then $𝔼 [𝑌_{𝑖}] = 0$ and $Var (𝑌_{𝑖}) = \frac{𝑉}{1} \approx \hat{𝑉}$ .

Then

\begin{aligned} err (𝑀 𝐶_{𝑁}) [𝑓, 𝑅] & = \frac{1}{𝑁} \sum_{𝑖} 𝑌_{𝑖} \\ ⟶ 𝑁 (0, \frac{𝑉}{𝑁}) by the central limit theorem \\ \approx 𝑁 (0, \frac{\hat{𝑉}}{𝑁}) by Lemma 3.2.6 . \end{aligned}

Corollary 3.2.9: The error of Monte Carlo integration is within 1 standard error about 68% of the time; within 2 standard errors 95% of the time; and within 3 standard error about 99.7% of the time.

Remark:

There is no guarantee on the error (you could get really unlucky)
Sampling even from strangely shaped regions is actually easy - just sample from a tight box around it, and reject any samples that are outside the region
This sampling technique could be quite time consuming in higher dimensions
Finding the area of the region could be quite difficult - but we can estimate it by considering the size of the sampling box and also how many samples were rejected. (Note that this is also a Monte Carlo integral! of the function ${\begin{cases} [1, 𝑥 \in 𝑅] \\ [0, 𝑥 \notin 𝑅] \end{cases}$ ). We therefore need to adjust the standard error.

Remark: Some improvements we could make:

Divide $𝑅$ into subregions, perhaps recursively (stratified sampling)
Sample some parts of $𝑅$ more densely, using the pdf of $𝑿_{𝒊}$ as $𝑔$ , then
$𝑀 𝐶_{𝑁} [𝑓, 𝑅] = \frac{1}{𝑁} \sum_{𝑖} \frac{𝑓 (𝑿_{𝒊})}{𝑔 (𝑿_{𝒊})}$
(importance sampling)
Subtract an easy integral (a control variate) to make the Monte Carlo integral smaller, so that the error is smaller

These don’t change the asymptotic error, but can give a better constant term.

4. Accuracy

Definition 4.1: Approximating $𝑢$ by $\tilde{𝑢}$ ,

$\tilde{𝑢} - 𝑢$ is the error
$| \tilde{𝑢} - 𝑢 |$ is the absolute error
$\frac{| \tilde{𝑢} - 𝑢 |}{| 𝑢 |}$ is the relative error.

For vector $𝒖$ , replace $| * |$ with the Euclidean norm $‖ * ‖$ .

Definition 4.2: If $𝑓 (𝑥)$ is an ideal function of $𝑥$ , and $\hat{𝑓} (𝑥)$ is an approximate implementation of $𝑓$ , then

$\hat{𝑓} (𝑥) - 𝑓 (𝑥)$ is the forward error
$\tilde{𝑥} - 𝑥$ is the backward error, where $\hat{𝑓} (𝑥) = 𝑓 (\tilde{𝑥})$ .

Definition 4.3: A floating-point number is a binary rational approximation to real numbers:

\tilde{𝑥} = \pm (1 + \frac{𝑚}{2^{𝑝}}) \cdot 2^{𝑒}

where $𝑚$ is the mantissa, $𝑝$ is the fixed precision (of the mantissa), and $𝑒$ is the exponent.

The IEEE 754 standard specifies:

single precision floats: 23 bit mantissa precision, 8 bit exponent, 1 sign bit
double precision: 52 bit mantissa precision, 11 bit exponent, 1 sign bit

Definition 4.4: Machine epsilon is the $𝜀$ such that every $𝑥 \in ℝ$ is representable with a relative error of no more than $𝜀$ .

Note that $𝜀$ is the relative error when we try to represent

𝑥 = 1 . \overset{𝑝 times}{\overset{⏞}{000 \dots 0}} 1_{2} \cdot 2^{𝑒}

for any representable exponent $𝑒$ .

Then

\begin{aligned} 𝜀 & = \frac{| 1 \cdot 2^{𝑒} - 1 . \overset{𝑝 times}{\overset{⏞}{000 \dots 0}} 1_{2} \cdot 2^{𝑒} |}{| 1.000 \dots 1 \cdot 2^{𝑒} |} \\ = \frac{| 2^{- (𝑝 + 1)} |}{| 1.000 \dots 1 |} \\ \approx 2^{- (𝑝 + 1)} . \end{aligned}

Machine epsilon is therefore $2^{- 24}$ for floats, and $2^{- 53}$ for doubles.

Some people define $𝜀$ to be the smallest $𝜀$ such that $1 + 𝜀$ is representable. This is double the definition used here.

Remark: Every

𝑥 \in ℝ

can be represented as a floating point number with relative error no more than

𝜀

Remark: If the result of a computation is orders of magnitude smaller than the inputs, we can get “catastrophic cancellation” leading to very large relative errors. E.g. subtraction of close numbers can round really badly for large inputs.

Definition 4.5: Truncation error occurs when approximating an infinite process by a finite process.

Roundoff error is due to floating-point storage and computation.

4.1. Iterative methods

Definition 4.1.1: Starting with an initial estimate $𝑥_{0}$ approximating unknown $𝑥^{*}$ , we repeatedly refine it with $𝑥_{𝑛 + 1} = 𝑓 (𝑥_{𝑛})$ , stopping when some criterion is met, e.g. $| 𝑥_{𝑛 + 1} - 𝑥_{𝑛} |$ < delta.

The error at step $𝑛$ is $𝜀_{𝑛} = 𝑥_{𝑛} - 𝑥^{*}$ .

If $𝜀_{𝑛} \to 0$ as $𝑛 \to \infty$ :

$𝑥_{𝑛}$ has linear convergence if for some $0 < 𝑎 < 1$ ,
$\frac{| 𝜀_{𝑛 + 1} |}{| 𝜀_{𝑛} |} \to 𝑎$
$𝑥_{𝑛}$ has sublinear convergence if
$| \frac{𝜀_{𝑛 + 1}}{| 𝜀 (𝑛) |} | \to 1$
- in particular $𝑥_{𝑛}$ has logarithmic convergence if
  $\frac{| 𝜀_{𝑛 + 2} - 𝜀_{𝑛 + 1} |}{| 𝜀_{𝑛 + 1} - 𝜀_{𝑛} |} \to 1$
$𝑥_{𝑛}$ converges superlinearly if
$\frac{| 𝜀_{𝑛 + 1} |}{| 𝜀_{𝑛} |} \to 0$
- in particular $𝑥_{𝑛}$ has order- $𝑞$ convergence for $𝑞 \geq 1$ if
  $\frac{| 𝜀_{𝑛 + 1} |}{{| 𝑛 |}^{𝑞}} \to 𝑎$
  for some $𝑎 > 0$

Remark: Linear convergence is called that not because the error is

𝑂 (𝑛)

, but because the amount of work you need to do to get an extra decimal place is constant, i.e. on a log plot we get a straight line.

Example:

$𝜀_{𝑛} = 2^{- 𝑛}$ converges linearly
$𝜀_{𝑛} = \frac{1}{𝑛^{𝑘}}$ converges logarithmically (a straight line on a log-log graph)
$𝜀_{𝑛} = 2^{- 2^{𝑛}}$ converges superlinearly, and in fact converges quadratically. On a log plot, this is an inverted parabola.

Remark: This can all be extended to vectors easily, by replacing

| * |

with

‖ * ‖

Remark: Before finding how quickly an iterative method converges, first determine if it converges at all: assume that it does converge, and show it converges to the right thing, then find an expression for

𝜀_{𝑛 + 1}

in terms of

𝜀_{𝑛}

, then in terms of

𝜀_{0}

with induction, and show that that tends to 0.

Lemma 4.1.2: If $𝑥_{𝑛}$ converges with order $𝑞$ , then $𝑥_{2 𝑛}$ converges with order $𝑞^{2}$ .

Proof: If $\frac{| 𝜀_{𝑛 + 1} |}{{| 𝜀_{𝑛} |}^{𝑞}} \to 𝑎$ for $𝑎 > 0$ , then

\frac{| 𝜀_{2 (𝑛 + 1)} |}{{| 𝜀_{2 𝑛} |}^{𝑞^{2}}} = \frac{| 𝜀_{2 𝑛 + 2} |}{{| 𝜀_{2 𝑛 + 1} |}^{𝑞}} {(\frac{| 𝜀_{2 𝑛 + 1} |}{{| 𝜀_{2 𝑛} |}^{𝑞}})}^{𝑞} \to 𝑎^{𝑞 + 1},

with

𝑎^{𝑞 + 1} > 0

5. Algorithms for root-finding

5.1. 1-dimensional root-finding

Definition 5.1.1: A root of

𝑓 : ℝ \to ℝ

𝑥^{*}

with

𝑓 (𝑥^{*}) = 0

Definition 5.1.2: A bracket is an interval $(𝑎_{0}, 𝑏_{0})$ with

𝑓 (𝑎_{0}) 𝑓 (𝑏_{0}) < 0 .

Theorem 5.1.3: If

𝑓

is continuous, every bracket contains at least one root.

Remark: This does not hold in higher dimensions.

Remark: We can effectively use binary search to find roots using brackets. However, we might not find an exact zero (because floating point), so we terminate when the interval is small enough, guess $𝑥^{*}$ is the midpoint of the interval, and can bound $𝑥^{*}$ by the bracket.

This is called interval bisection.

Iterating, we get a sequence $𝑥_{0}, 𝑥_{1}, \dots$ with

| 𝜀_{𝑛} | < \frac{| 𝑏_{0} - 𝑎_{0} |}{2^{𝑛 + 1}}

so this method converges linearly.

Definition 5.1.4: Newton’s method in 1 dimension approximates $𝑓$ by its first-order Taylor polynomial.

\hat{𝑓} (𝑥) = 𝑓 (𝑥_{0}) + (𝑥 - 𝑥_{0}) \frac{d 𝑓}{d 𝑥} (𝑥_{0}) .

Take a guess $𝑥_{0}$ , then find the root of the first-order Taylor polynomial at $𝑥_{0}$ ; this is $𝑥_{1}$ , i.e.

\begin{matrix} 0 = 𝑓 (𝑥_{𝑛}) + (𝑥_{𝑛 + 1} - 𝑥_{𝑛}) \frac{d 𝑓}{d 𝑥} (𝑥_{𝑛}) \\ ⟹ 𝑥_{𝑛 + 1} = 𝑥_{𝑛} - \frac{𝑓 (𝑥_{𝑛})}{\frac{d 𝑓}{d 𝑥} (𝑥_{𝑛})} . \end{matrix}

This converges quadratically for some $𝑥_{0}$ . We can also get stuck at a fixed point, or hit a turning point and divide by zero, or diverge.

Remark: Some limitations of Newton’s method:

you need to be able to evaluate the derivative
fails if $\frac{d 𝑓}{d 𝑥} (𝑥_{𝑛}) = 0$ (throw an error)
iteration can get stuck in a loop
iteration can diverge

Theorem 5.1.5: The forward error estimate for Newton’s method is

| 𝑥^{*} - 𝑥_{𝑛} | \approx | 𝑥_{𝑛 + 1} - 𝑥_{𝑛} | .

Proof:

\begin{aligned} \hat{𝑓} (𝑥^{*}) & = 𝑓 (𝑥_{𝑛}) + (𝑥^{*} - 𝑥_{𝑛}) \frac{d 𝑓}{d 𝑥} (𝑥_{𝑛}) \\ \approx 0 for root 𝑥^{*} . \\ ∴ | 𝑥_{𝑛} - 𝑥^{*} | & \approx | \frac{𝑓 (𝑥_{𝑛})}{\frac{d 𝑓}{d 𝑥} (𝑥_{𝑛})} | \\ = | 𝑥_{𝑛 + 1} - 𝑥_{𝑛} | by the definition of 𝑥_{𝑛 + 1} . \end{aligned}

Remark: So we might choose to terminate when successive iterations are sufficiently close, by some parameter $𝑡 𝑜 𝑙$ , often $2 𝜀$ or $4 𝜀$ , preferably by some relative error.

We can also look at backward error, $| 𝑓 (𝑥_{𝑛}) | < 𝑡 𝑜 𝑙$ .

There’s no one best way of doing a stopping condition, you should think about it.

Remark: You also must stop at a sensible upper bound of steps so you don’t get stuck in an infinite loop (ideally with a warning that convergence didn’t happen).

You must also stop if an iteration step is poorly defined.

Lemma 5.1.6:

Let $𝐼_{𝑐} = (𝑥^{*} - 𝑐, 𝑥^{*} + 𝑐)$ for some $𝑐 \in ℝ$ .

Let

𝐴 (𝑐) = \frac{{max}_{𝛽 \in 𝐼_{𝑐}} | \frac{d^{2} 𝑓}{d 𝑥^{2}} (𝛽) |}{{min}_{𝛼 \in 𝐼_{𝑐}} | \frac{d 𝑓}{d 𝑥} (𝛼) |} .

Then if $𝑥_{𝑛} \in 𝐼_{𝑐}$ , then $| 𝜀_{𝑛 + 1} | \leq \frac{𝐴 (𝑐)}{2} 𝜀_{𝑛}^{2}$ .

Proof:

Suppose $𝑥_{𝑛} \in 𝐼_{𝑐}$ .

\begin{aligned} | 𝜀_{𝑛 + 1} | & = | 𝑥_{𝑛 + 1} - 𝑥^{*} | \\ = | 𝑥_{𝑛} - \frac{𝑓 (𝑥_{𝑛})}{𝑓^{'} (𝑥_{𝑛})} - 𝑥^{*} | \\ = | \frac{(𝑥_{𝑛} - 𝑥^{*}) 𝑓^{'} (𝑥_{𝑛}) - 𝑓 (𝑥_{𝑛})}{𝑓^{'} (𝑥_{𝑛})} | \\ = \frac{| (𝑓 (𝑥^{*}) - (𝑓 (𝑥_{𝑛}) + (𝑥^{*} - 𝑥_{𝑛}) 𝑓^{'} (𝑥_{𝑛}))) |}{| (𝑓^{'} (𝑥_{𝑛})) |} because 𝑓 (𝑥^{*}) = 0 \\ = \frac{| \frac{1}{2} {(𝑥_{𝑛} - 𝑥^{*})}^{2} 𝑓^{″} (𝜉) |}{| 𝑓^{'} (𝑥_{𝑛}) |} for some 𝜉 \in (𝑥^{*}, 𝑥_{𝑛}), by Taylor's theorem \\ = \frac{1}{2} 𝜀_{𝑛}^{2} \frac{| 𝑓^{″} (𝜉) |}{| 𝑓^{'} (𝑥_{𝑛}) |} \\ \leq \frac{𝐴 (𝑐)}{2} 𝜀_{𝑛}^{2} . \end{aligned}

Lemma 5.1.7:

For $𝐴 (𝑐), 𝐼_{𝑐}$ as before, if $𝑥_{0} \in 𝐼_{𝑐}$ and $\frac{𝑐 𝐴 (𝑐)}{2} < 1$ then $𝑥_{𝑛}$ converges to $𝑥^{*}$ at least quadratically.

Proof:

Suppose $\frac{𝑐 𝐴 (𝑐)}{2} < 1$ and that $𝑥_{𝑛} \in 𝐼_{𝑐}$ as before. Note that $| 𝜀_{𝑛} | < 𝑐$ from the definition of $𝐼$ .

\begin{aligned} | 𝜀_{𝑛 + 1} | & \leq \frac{𝐴 (𝑐)}{2} 𝜀_{𝑛}^{2} from Lemma 5.1.6 \\ = | 𝜀_{𝑛} | \frac{| 𝜀_{𝑛} | 𝐴 (𝑐)}{2} \\ \leq | 𝜀_{𝑛} | because | 𝜀_{𝑛} | \frac{𝐴 (𝑐)}{2} < 𝑐 \frac{𝐴 (𝑐)}{2} < 1 by supposition \\ ∴ 𝑥_{𝑛 + 1} & \in 𝐼 . \end{aligned}

Then letting $𝜌 = \frac{| 𝜀_{𝑛} | 𝐴 (𝑐)}{2}$ , by induction, $𝑥_{𝑛} \in 𝐼$ for all $𝑛$ , and $| 𝜀_{𝑛} | \leq 𝜌^{𝑛} | 𝑥_{0} |$ . $𝜌 < 1$ so $| 𝜀_{𝑛} | \to 0$ as $𝑛 \to \infty$ .

Therefore $𝑥_{𝑛} \to 𝑥^{*}$ .

Moreover,

\frac{| 𝜀_{𝑛 + 1} |}{{| 𝜀_{𝑛} |}^{2}} \leq \frac{𝐴 (𝑐)}{2},

so we have at least quadratic convergence (it might be better if it actually tends to 0; we don’t have equality).

Remark: The important thing here is that we have proven that the

𝑥_{𝑛}

s stay within

𝐼_{𝑐}

Remark: The condition on

𝑐

here,

\frac{𝑐 𝐴 (𝑐)}{2} < 1

, is a sufficient condition for

𝐼_{𝑐}

to be a basin of convergence.

Lemma 5.1.8: Let $𝐽_{𝑐} = (𝑥_{0} - 𝑐, 𝑥_{0} + 𝑐)$ , and let

𝐵 (𝑐) = \frac{{max}_{𝛽 \in 𝐽_{𝑐}} | \frac{d^{2} 𝑓}{d 𝑥^{2}} (𝛽) |}{{min}_{𝛼 \in 𝐽_{𝑐}} | \frac{d 𝑓}{d 𝑥} (𝛼) |} .

Then if $𝑥^{*} \in 𝐽_{𝑐}$ and $\frac{𝑐 𝐵 (𝑐)}{2} < 1$ , then $𝑥_{𝑛} \to 𝑥^{*}$ quadratically.

Proof:

Suppose that $𝑥^{*} \in 𝐽_{𝑐}$ and $\frac{𝑐 𝐵 (𝑐)}{2} < 1$ .

From Lemma 5.1.6,

| 𝜀_{𝑛 + 1} | \leq 𝜀_{𝑛}^{2} \frac{𝐴 (𝑐)}{2} < | 𝜀_{𝑛} |,

because the interval contained within $𝐽_{𝑐}$ with $𝑥^{*}$ at its centre, $𝐼_{𝑑}$ , must have $𝑑 \leq 𝑐$ , so $𝐴 (𝑑) \leq 𝐴 (𝑐)$ .

Note that although we cannot guarantee that $𝑥_{𝑛}$ stays within $𝐽_{𝑐}$ ,we can guarantee that it stays in $𝐽_{2 𝑐}$ , if $\frac{𝑐 𝐵 (𝑐)}{2} < 1$ , by the same reasoning as Lemma 5.1.7.

TODO make this a bit nicer from lecture notes

Remark: We might find such an

𝑥_{0}, 𝑐

by finding a suitable bracket.

Theorem 5.1.9: There exists an $𝑥_{0}, 𝑐$ such that $𝐽_{𝑐}$ is a basin of convergence if $\frac{d 𝑓}{d 𝑥} (𝑥^{*}) \neq 0$ .

Proof: Given a bracket containing

(𝑥^{*} - 𝑐, 𝑥^{*} + 𝑐)

, we can shrink the interval down such that

𝐴 (𝑐) \to \frac{𝑓^{″} (𝑥^{*})}{𝑓^{'} (𝑥^{*})}

. If

𝑓^{'} (𝑥^{*}) \neq 0

, then this is finite and we can find small

𝑐

such that

\frac{1}{2} 𝑐 𝐴 (𝑐) < 1

; the same is true for

𝐵 (𝑐)

Remark: The convergence result doesn’t work for double roots (it might still converge, but only linearly).

We can fix it by using the iteration $𝑥_{𝑛 + 1} = 𝑥_{𝑛} - 2 \frac{𝑓 (𝑥_{𝑛})}{𝑓^{'} (𝑥_{𝑛})}$ , which will converge quadratically to double roots.

Definition 5.1.10: The secant method uses a linear interpolation to approximate $𝑓$ .

𝑥_{𝑛 + 2} = 𝑥_{𝑛 + 1} - \frac{𝑓 (𝑥_{𝑛 + 1}) (𝑥_{𝑛 + 1} - 𝑥_{𝑛})}{𝑓 (𝑥_{𝑛 + 1}) - 𝑓 (𝑥_{𝑛})} .

Remark: This is a “quasi-Newton method”, where we use an estimate for the derivative.

Remark: This only needs one evaluation of

𝑓

at each iteration.

Remark: This can converge superlinearly, in fact order-

𝜙

in good cases. (proof: exercise)

Remark: In theory this can diverge or get stuck in a loop, but it is more stable than Newton’s method, because it relies on the previous 2 iterations.

TODO forward error (same as Newton)

Remark: We should be able to do two iterations of the secant method for no more than the cost of 1 iteration of Newton’s method, which would converge with order

\approx 2.6

. This is also nice because we just don’t want to evaluate the derivative of

𝑓

Remark: Some other 1-dimensional root finding methods:

Halley’s method - higher order Taylor approximation. But difficult because we might end up with no roots or many roots. But cubic convergence in good cases.
Muller’s method - higher order polynomial interpolation - same problems as Halley’s method, but no derivatives needed and order $\approx 1.84$ when it works
Inverse quadratic interpolation, i.e. find $𝑥 = 𝑎 𝑦^{2} + 𝑏 𝑦 + 𝑐$ , and the root is at $𝑐$ (always!). Order $\approx 1.84$ convergence in good cases.

But the best method is Brent’s method:

maintains a bracket
uses inverse quadratic interpolation unless iteration undefined or outside the bracket
uses the secant method as second choice
falls back to bisection if these do not sufficiently shrink the bracket (only happens at at most half of iterations)
derivative free
always converges, at order $\approx 1.84$ in good cases.

5.2. $𝑑$ -dimensional root finding

Definition 5.2.1: A root of $𝒇 : ℝ^{𝑑} \to ℝ^{𝑘}$ is $𝑿^{*}$ with $𝒇 (𝑿^{*}) = 𝟎$ .

If $𝑘 < 𝑑$ , the problem is underconstrained, so there are probably infinitely many solutions.

If $𝑘 > 𝑑$ , the problem is unconstrained, so there are probably no (exact) solutions.

We will consider $𝑘 = 𝑑$ .

Remark: There isn’t really such a thing as a bracket in higher dimensions, because although we might have roots of components, we might not find simultaneous roots of these.

Remark: If instead we take a hyperrectangle with a bracket for components on opposite sides, then we can guarantee a simultaneous root.

TODO diagram

But this is not very useful, because we can’t algorithmically prove tha a function is positive/negative all the way along a line.

Remark: For higher-dimensional iterative methods, the sensible termination conditions are the same as they were in one dimension, but with norms rather than absolute value, and a tolerance parameter should probably rely on

𝑑

Definition 5.2.2: Newton’s method for $𝑑$ -dimensional root-finding can be derived as follows.

Given $𝒇 : ℝ^{𝑑} \to ℝ^{𝑑}, 𝒇 (𝒙) = (\begin{matrix} 𝑓_{1} (𝒙) \\ \dots \\ 𝑓_{𝑑} (𝒙) \end{matrix})$ , we approximate each component by its first-order Taylor polynomial:

{\hat{𝑓}}_{𝑖} (𝒙) = 𝑓_{𝑖} (𝒙_{𝟎}) + {(\frac{d 𝑓_{𝑖}}{d 𝒙} (𝒙_{𝟎}))}^{⊤} (𝒙 - 𝒙_{𝟎}) .

Then we want to find roots of all of these approximations:

\begin{matrix} 0 = 𝑓_{𝑖} (𝒙_{𝟎}) + {(\frac{d 𝑓_{𝑖}}{d 𝒙} (𝒙_{𝟎}))}^{⊤} (𝒙 - 𝒙_{𝟎}) \\ ⟹ {(\frac{d 𝑓_{𝑖}}{d 𝒙} (𝒙_{𝟎}))}^{⊤} (𝒙 - 𝒙_{𝟎}) = - 𝑓_{𝑖} (𝒙_{𝟎}) . \end{matrix}

We can put all of these together to get

𝐉 (𝒇) (𝒙_{𝟎}) (𝒙 - 𝒙_{𝟎}) = - 𝒇 (𝒙_{𝟎}) .

So we get an iteration

𝒙_{𝑛 + 1} = 𝒙_{𝑛} - 𝐉 (𝒇) {(𝒙_{𝑛})}^{- 1} 𝒇 (𝒙_{𝑛}) .

Each iteration, this requires 1 evaluation of $𝑓$ , $𝑑^{2}$ derivative, and an inversion, which is $𝑂 (𝑑^{3})$ .

However note that actually just solving the first equation with the Jacobian, rather than inverting it, is a constant factor quicker.

Remark: This is not the curse of dimensionality; it only blows up polynomially with

𝑑

, not exponentially.

Remark: We still have existence of a basin of quadratic convergence, if:

$𝒇, 𝐉 (𝒇), \frac{𝜕^{2} 𝑓_{𝑖}}{𝜕 𝑥_{𝑗} 𝜕 𝑥_{𝑘}}$ are continuous
$𝐉 (𝒇)$ is nonsingular at $𝒙^{*}$

Remark: In order to globalise Newton’s method so that it works everywhere, we should enforce that if $‖ 𝒇 (𝒙_{𝑛 + 1}) ‖ < ‖ 𝒇 (𝒙_{𝑛}) ‖$ , we should backtrack by setting $𝒙_{𝒏 + 𝟏} = 𝒙_{𝒏} + 𝜆 Δ 𝒙$ , where $Δ 𝒙$ is the $(𝒙 - 𝒙_{𝟎})$ that is solved at each step, for some $𝜆 \in (0, 1)$ .

Suitably small $𝜆$ guarantees linear convergence everywhere, and this is called a damped Newton’s method.

Definition 5.2.3: In $𝑑$ dimensions, the equivalent to the secant method is Broyden’s method. On each step we use an approximation for the Jacobian, ${\hat{𝐉}}_{𝑛}$ .

We would like for $𝑓 (𝒙_{𝒏}) = 𝑓 (𝒙_{𝒏 - 𝟏}) + {\hat{𝐉}}_{𝑛} (𝒙_{𝑛} - 𝒙_{𝒏 + 𝟏})$ , but we cannot solve for ${\hat{𝐉}}_{𝑛}$ because it is an underdetermined system.

We can notice that our desired relation doesn’t tell us anything about the relation between ${\hat{𝐉}}_{𝑛}$ and vector $𝑧$ that is orthogonal to $(𝒙_{𝒏} - 𝒙_{𝒏 - 𝟏})$ ; we could therefore choose that consecutive ${\hat{𝐉}}_{𝑛}$ s should behave in the same manner for such vectors, i.e.

{\hat{𝐉}}_{𝑛} 𝒛 = {\hat{𝐉}}_{𝑛 - 1} 𝒛, for 𝒛^{⊤} (𝒙_{𝒏} - 𝒙_{𝒏 - 𝟏}) = 0 .

We now have enough constraints to find ${\hat{𝐉}}_{𝑛}$ , provided we have an initial ${\hat{𝐉}}_{0}$ , which we could either calculate as $𝐉 (𝑓) (𝒙_{𝟎})$ , or set to $- 𝐈$ , or a multiple of.

Equivalently, we could aim to minimise the norm of the difference between consecutive Jacobian estimates, which makes sense because close to the root, they should be similar.

Then, writing $𝚫 𝒙 = 𝒙_{𝒏} - 𝒙_{𝒏 - 𝟏}$ ,

{\hat{𝐉}}_{𝑛} = {\hat{𝐉}}_{𝑛 - 1} + \frac{(𝑓 (𝒙_{𝒏}) - 𝑓 (𝒙_{𝒏 - 𝟏})) - {\hat{𝐉}}_{𝑛} 𝚫 𝒙}{{‖ 𝚫 𝒙 ‖}^{2}} {𝚫 𝒙}^{⊤} .

This step has complexity $𝑂 (𝑑^{2})$ because the most we do is matrix-vector multiplication.

Then the iteration is

𝒙_{𝒏 + 𝟏} = 𝒙_{𝒏} + 𝚫 𝒙, where {\hat{𝐉}}_{𝑛} 𝚫 𝒙 = - 𝑓 (𝒙_{𝒏}) .

This step still requires $𝑂 (𝑑^{3})$ to solve the system, so overall each iteration is $𝑂 (𝑑^{3})$ .

The convergence is order $𝑞 > 1$ in good cases.

Theorem 5.2.4:

{(𝐀 + 𝒖 𝒗^{⊤})}^{- 1} =

TODO

Proof: Exercise

Remark: This means that we can start from ${\hat{𝐉}}_{0}^{- 1}$ and never need to calculate any inverses.

This gives $𝑂 (𝑑^{2})$ for Broyden’s method.

6. Algorithms for optimisation

6.1. 1-dimensional minimisation

Remark: If we can show that a function is convex, we can use a root finding algorithm to find the minimum. We can find a bracket for the minimum by finding a bracket for a root of the derivative (but for a minimum it;s got to be negative on the left and positive on the right).

Definition 6.1.1: We can find a bracket for a minimum without evaluating the derivative: if there is $𝑎 < 𝑧 < 𝑏$ s.t. $𝑓 (𝑧) < 𝑓 (𝑎), 𝑓 (𝑧) < 𝑓 (𝑏)$ , then $(𝑎, 𝑏)$ is a bracket containing the minimum.

We can make an algorithm for reducing the bracket, by keeping track of the middle value and testing at a new point $𝑧^{'}$ , and discarding one of the endpoints. This is called golden section search.

Let’s ignore the case of a tie between the middle values, and suppose w.l.o.g. that we choose $𝑧^{'} \in (𝑧, 𝑐)$ (we want to choose from the larger side of the bracket so that we remove us much as possible, and the other case is symmetric).

The placement of $𝑧^{'}$ is quite important: we want ideally the same amount to be chopped off whether $𝑧^{'} > 𝑧$ or not; therefore we say that $𝑧^{'} - 𝑎 = 𝑐 - 𝑧$ .

We would also like the bracket to shrink by a constant factor each time, so we want $𝑐 - 𝑧^{'} = 𝜙 (𝑐 - 𝑧)$ , and also $𝑐 - 𝑧 = 𝜙 (𝑐 - 𝑎)$ .

Then

\begin{aligned} 𝑐 - 𝑎 & = 𝑐 - 𝑧^{'} + 𝑧^{'} - 𝑎 \\ ⟹ \frac{1}{𝜙^{2}} (𝑐 - 𝑧^{'}) & = (𝑐 - 𝑧^{'}) + (𝑐 - 𝑧) = (𝑐 - 𝑧) + \frac{1}{𝜑} (𝑐 - 𝑧^{'}) . \end{aligned}

Therefore $\frac{1}{𝜙^{2}} = 1 + \frac{1}{𝜙} ⟹ 𝜙^{2} + 𝜙 - 1 = 0$ ; that is, $𝜙$ is the golden ratio, hence the name golden section search.

Proposition 6.1.2: This converges linearly.

Proof: TODO

Remark: Our termination condition should be

| 𝑐 - 𝑎 | \leq 𝑡 𝑜 𝑙 | 𝑧 |

, where

𝑡 𝑜 𝑙

should only be about

𝑂 (\sqrt{𝜀})

Definition 6.1.3: Successive parabolic interpolation carries out golden section search but by choosing $𝑧^{'}$ to be the minimum of the parabola passing through $𝑎, 𝑧, 𝑏$ .

This is order $\approx 1.32$ , but it can go wrong in many ways (exercise).

Remark: Finding a bracket isn’t necessarily easy, but if the function is convex, we can use a geometrically expanding bracket to find a bracket. Even better, we can keep the golden ratio, with the geometrically expansion.

Definition 6.1.4: Brent’s method for minimisation:

combines successive parabolic interpolation with golden section search, keeping track of 6 points
is derivative-free and superlinear in good cases.

6.2. $𝑑$ -dimensional minimisation

Given a position $𝒙_{𝒏}$ , choose direction $𝒅_{𝒏}$ and step length $𝛼_{𝑛}$ , and set $𝒙_{𝒏 + 𝟏} = 𝒙_{𝒏} + 𝛼_{𝑛} 𝒅_{𝒏}$ .

Define $𝒈_{𝒏} = \frac{d 𝑓}{d 𝒙} (𝒙_{𝒏})$ , with $𝒈$ standing for gradient.

We would like $𝒅_{𝒏}$ to be a descent direction:

\begin{aligned} \frac{d 𝑓 (𝒙_{𝒏} + 𝛼_{𝑛} 𝒅_{𝒏})}{d 𝛼} < 0 \\ ⟹ & {(\frac{d 𝑓}{d 𝒙} \circ (𝛼 \mapsto 𝒙_{𝒏} + 𝛼 𝒅_{𝒏}))}^{⊤} 𝐉_{𝛼} (𝒙_{𝒏} + 𝛼 𝒅_{𝒏}) < 0 \\ ⟹ & \frac{d 𝑓}{d 𝒙} {(𝒙_{𝒏} + 𝛼 𝒅_{𝒏})}^{⊤} 𝒅_{𝒏} < 0 \\ ⟹ & {𝒈_{𝒏}}^{⊤} 𝒅_{𝒏} < 0 . \end{aligned}

We’d also like consecutive directions shouldn’t veer wildly, undoing work of previous steps; they should perhaps be roughly orthogonal, to explore the space efficiently.

$𝛼_{𝑛}$ should loosely minimise $𝑓 (𝒙_{𝒏} + 𝛼_{𝑛} 𝒅_{𝒏})$ .

The “infinite travel” condition should also be met, to avoid the process halting completely:

\sum_{𝑛 = 1}^{\infty} 𝛼_{𝑛} = \infty .

Definition 6.2.1: Coordinate gradient descent chooses alternating basis vectors for the direction (with

𝛼_{𝑛}

then chosen to be negative if it is in the wrong direction). It works, but it’s not very good.

Definition 6.2.2: Gradient descent finds the steepest downhill direction:

\begin{matrix} min {𝒅_{𝒏}}^{⊤} 𝒈_{𝒏} s.t. 1 - ‖ 𝒅_{𝑛} ‖ \geq 0 . \\ Λ (𝒅_{𝒏}, 𝜇) = {𝒅_{𝒏}}^{⊤} 𝒈_{𝒏} - 𝝁 (1 - {𝒅_{𝒏}}^{⊤} 𝒅_{𝒏}) \\ \frac{d Λ}{d 𝒅_{𝒏}} = 𝒈_{𝒏} + 2 𝜇 𝒅_{𝒏} = 𝟎 \\ ⟹ 𝒅 = - \frac{𝒈_{𝒏}}{2 𝜇} \end{matrix}

So we choose $𝒅 = - 𝒈_{𝒏}$ , as the scalar multiplication can be chosen by $𝛼_{𝑛}$ .

That is, we just go in the direction of the negative gradient; this makes sense.

To find a good step length, overshoot and backtrack (typically halving) until the step length is acceptable.

Our acceptable condition for when we can stop backtracking is

𝑓 (𝒙_{𝒏} + 𝛼_{𝑛} 𝒅_{𝒏}) < 𝑓 (𝒙_{𝒏}) + 𝜎 𝛼_{𝑛} {𝒈_{𝒏}}^{⊤} 𝒅_{𝒏} .

This is called the Armijo rule. $𝜎$ is typically $0.1$ to $0.0001$ . It says that we want the new value of the function to be below a line that is somewhere in between the 0th order and 1st order approximation to f at $𝑥_{𝑛}$ .

It can be proven that this always works if $𝜎 < 1$ .

We could take an initial guess for $𝛼_{𝑛}$ to be

𝛼_{𝑛 - 1} \frac{{𝒈_{𝒏 - 𝟏}}^{⊤} 𝒅_{𝒏 - 𝟏}}{{𝒈_{𝒏}}^{⊤} 𝒅_{𝒏}} .

i.e. we suppose that an approximation for the amount of descent in this step will be the same as in the previous step. It is better to overestimate the initial $𝛼_{𝑛}$ than to underestimate, as backtracking is cheap.

We still need a large initial step length for step 0. We could even forward-track until we get something that doesn’t meet the Armijo rule.

TODO overall iteration

TODO termination condition.

TODO complexity

TODO convergence

Definition 6.2.3: Conjugate gradient descent uses conjugate directions that satisfy $𝒅_{𝒏} 𝐇 𝒅_{𝒏} = 0$ .

𝒅_{𝒏} = - 𝒈_{𝒏} + 𝒅_{𝒏 - 𝟏} (

TODO)TODO

TODO justification, what does this mean?

Remark: Conjugate gradient descent isn’t actually on the syllabus but it’s useful and quite simple (only one line of code changed compared to gradient descent).

Remark: We can also use Newton’s method:

Approximate $𝑓$ by its second-order Taylor polynomial

\hat{𝑓} (𝒙) = 𝑓 (𝒙_{𝟎}) + {(𝒙 - 𝒙_{𝟎})}^{⊤} 𝒈_{𝟎} + \frac{1}{2} (𝒙 - 𝒙_{𝟎}) 𝐇 (𝑓) (𝒙_{𝟎}) (𝒙 - 𝒙_{𝟎}) .

Differentiating this gives

\frac{d \hat{𝑓}}{d 𝒙} = 𝟎 + 𝒈_{𝒏} - 𝟎 + 𝐇_{0} (𝒙 - 𝒙_{𝟎}) .

where $𝐇_{𝑛} = 𝐇 (𝑓) (𝒙_{𝒏})$ .

Therefore $𝒙_{𝟏} - 𝒙_{𝟎} = - 𝐇_{0}^{- 1} 𝒈_{𝟎}$ . But remember that computing matrix inverses is slow, so better to say solve $𝐇_{0} 𝒅_{𝒏} = - 𝒈_{𝒏}$ for $𝒅_{𝒏}$ .

This is like gradient descent, but with $𝒅_{𝒏} = - {(𝐇 (𝑓) (𝒙_{𝒏}))}^{- 1} 𝒈_{𝒏}$ (we call this the Newton direction) and $𝛼_{𝑛} = 1$ .

Note that this has problem-independent step length, which is nice.

However, the requirement for computing the Hessian at every iteration makes each step a lot more expensive than gradient descent.

Remark: Newton’s method doesn’t always give a descent direction.

𝑔_{𝑛}^{⊤} 𝑑_{𝑛} = - 𝑔_{𝑛}^{⊤} 𝐻_{𝑛}^{- 1} 𝑔_{𝑛} < 0 if 𝐻_{𝑛}^{- 1} pos def ⟸ 𝐻_{𝑛} pos def .

Therefore Newton’s method can converge to maxima or saddle points equally as much as minima.

TODO compare the two Newton’s methods (they’re the same)

Remark: We could use a quasi-Newton method whereby we use an approximation of the Hessian, e.g. by using Broyden’s method on the derivative. However, this is bad because it doesn’t necessarily give a symmetric approximation for the Hessian.

We want ${\hat{𝐇}}_{𝑛}$ :

to satisfy secant equation TODO
to be symmetric
to be close to ${\hat{𝐇}}_{𝑛 - 1}$
to be positive definite.

Positive definiteness requires the curvature condition $Δ 𝒙^{⊤} Δ 𝒈 > 0$ , which we can’t always guarantee.

Definition 6.2.4: BFGS update is TODO (not in spec).

If the curvature condition fails, reset ${\hat{𝐇}}_{𝑛}$ .

Best to start with ${\hat{𝐇}}_{0} = 𝐈$ , and reset to the identity if the curvature condition fails. When the approximate Hessian is the identity, it is just gradient descent; so it does gradient descent until it gets to a convex section of the function, at which point it starts approximating the Hessian and goes a lot quicker.

TODO advantages

TODO disadvantages

Remark: There is a way to only update the inverse Hessians, like we did for Broyden’s method. It’s a bit more complicated because this uses a rank-2 update rather than a rank-1 update.

Remark: We can make a rank-

𝑚

approximation of the Hessians and get

𝑂 (𝑚 𝑑)

memory. We call this L-BFGS.

HT2026 Continuous Maths notes

1. Derivatives and Taylor’s Theorem

1.1. Continuity and differentiation

1.2. Taylor’s theorem for univariate functions

1.3. Partial and vector derivatives

1.4. Vector fields and the Jacobian

1.5. Taylor’s theorem for multivariate functions

2. Optimisation

2.1. Optimisation in 1 dimension

2.2. Positive/negative (semi)definiteness

2.3. Unconstrained optimisation over ℝ𝑛

2.4. Convexity

2.5. Optimisation tricks

2.6. Optimisation over ℝ𝑛 with equality constraints

2.7. Inequality contrained optimisation (over ℝ𝑛)

3. Algorithms for numerical integration

3.1. 1-dimensional integration

3.2. Integration in 𝑑 dimensions

4. Accuracy

4.1. Iterative methods

5. Algorithms for root-finding

5.1. 1-dimensional root-finding

5.2. 𝑑-dimensional root finding

6. Algorithms for optimisation

6.1. 1-dimensional minimisation

6.2. 𝑑-dimensional minimisation

2.3. Unconstrained optimisation over $ℝ^{𝑛}$

2.6. Optimisation over $ℝ^{𝑛}$ with equality constraints

2.7. Inequality contrained optimisation (over $ℝ^{𝑛}$ )

3.2. Integration in $𝑑$ dimensions

5.2. $𝑑$ -dimensional root finding

6.2. $𝑑$ -dimensional minimisation