Damped Newton’s Method with Self-Concordant Function

Ziwen Wang

February 16, 2026

1 Theoretical discussion

In this writing we will review how damped newton method works and provide an implementation of it. Throughout, we will assume that the given function $f : ℝ^{m} \to ℝ$ is strictly convex and twice differentiable. Newton’s method can be thought of finding the minimizer of Taylor’s second order approximation,

f (x) \approx f (x_{k}) + ⟨ x - x_{k}, \nabla f (x_{k}) ⟩ + \frac{1}{2} (x - x_{k}) \nabla^{2} f (x_{k}) (x - x_{k}) .

The exact minimizer assuming strict convexity is unique and has the following form,

{arg}_{x \in ℝ^{m}} {f (x_{k}) + ⟨ x - x_{k}, \nabla f (x_{k}) ⟩ + \frac{1}{2} (x - x_{k}) \nabla^{2} f (x_{k}) (x - x_{k})} = x_{k} - {(\nabla^{2} f (x_{k}))}^{- 1} \nabla f (x) .

Let us define the decrement direction at iteration $k$ as $d_{k} = - {(\nabla^{2} f (x_{k}))}^{- 1} \nabla f (x)$ , then the update step can be rewritten as,

x_{k + 1} = x_{k} + α_{k} d_{k}, α_{k} = 1 .

With step size $α_{k} = 1$ , this update rule is referred as the pure newton step, while this seems fine but it only converges (in fact quadratically) if the initial starting point $x_{0}$ is close enough to minimizer and with the assumption that the function $f$ is strongly convex and Lipschtiz continuous.

Lemma 1 (Newton’s method is affine invariant). Given an invertible matrix $A \in ℝ^{m \times m}$ and a convex function $f : ℝ^{m} \to ℝ$ . Let $x_{i}, {\tilde{x}}_{i}$ be newton iterates of $f (x), \tilde{f} (x) : = f (Ax)$ with starting point $x_{0}, {\tilde{x}}_{0} : = A^{- 1} x_{0}$ respectively. Then ${\tilde{x}}_{i} = A^{- 1} x_{i}$ .

Proof. Notice that $\nabla_{\tilde{x}} \tilde{f} (\tilde{x}) = A^{T} \nabla_{x} f (Ax)$ and $\nabla_{\tilde{x}}^{2} \tilde{f} (Ax) = A^{T} \nabla_{x}^{2} f (Ax) A$ and by induction we have

\begin{array}{l} {\tilde{x}}_{i + 1} & = {\tilde{x}}_{i} - {(\nabla^{2} f ({\tilde{x}}_{i}))}^{- 1} \nabla f ({\tilde{x}}_{i}) \\ = {\tilde{x}}_{i} - A^{- 1} \nabla^{2} f (x_{i}) A^{- T} A^{T} \nabla f (x_{i}) \\ = {\tilde{x}}_{i} - A^{- 1} \nabla^{2} f (x_{i}) \nabla f (x_{i}) \\ = A^{- 1} x_{i} - A^{- 1} \nabla^{2} f (x_{i}) \nabla f (x_{i}) \\ = A^{- 1} x_{i + 1} . \end{array}

□

Note at optimality we have $\nabla f (x) = 0$ , so one way to measure convergence is by bounding upper bounding $∥ f (x) ∥$ , but notice that $∥ f (x) ∥$ is not affine invariant i.e. $∥ f (x) ∥ \neq ∥ f (Ax) ∥$ so instead we have the following measure,

λ_{f} (x) : = \sqrt{{(\nabla f (x))}^{T} {(\nabla^{2} f (x))}^{- 1} \nabla f (x)} = : ∥ \nabla f (x) ∥_{{(\nabla^{2} f (x))}^{- 1}}

which becomes affine invariant and can be interpreted as an induced norm by the PD matrix ${(\nabla^{2} f (x))}^{- 1}$ . So now instead of bounding on $∥ f (x) ∥$ we can try to bound $λ_{f} (x)$ .

Under the assumption that $f$ is self concordant with constant 1 one can show that $λ_{f} (x_{k + 1}) \leq {\frac{λ_{f} (x_{k})}{1 - λ_{f} (x_{k})}}^{2}$ (in order for it to decrease we need $λ_{f} (x_{k}) < 1$ ), but this still assumes a good starting point (which is the condition $λ_{f} (x_{k}) < 1$ ). It turns out if we change the step size then we are able to achieve global convergence but before we do so, we shall define self-concordant functions more formally.

Definition 2 (Self-concordant function). A function $f$ is called self-concordant if there exists a constant $M_{f} > 0$ such that

| D^{3} f (x) [u, u, u] | \leq 2 M_{f} ∥ u ∥_{\nabla^{2} f (x)}^{3}

holds for all $x \in domf \subseteq ℝ^{m}$ and $u \in ℝ^{m}$ .

We now give three different stepping sizes and discuss their convergence properties, their advantage and disadvantages.

1.: (Standard Newton’s method) $α_{k} = 1$ ,
2.: (Damped Newton’s method) $α_{k} = \frac{1}{1 + M_{f} λ_{k}}$ ,
3.: (Intermediate Newton’s method) $α_{k} = \frac{1}{1 + \frac{{(M_{f} λ_{k})}^{2}}{1 + M_{f} λ_{k}}}$ , where $λ_{k} : = λ_{f} (x_{k})$ .

From above it is not hard to see that the following inequality of the steps sizes holds,

(1 .) \geq (3 .) \geq (2 .)

Method $(1 .)$ uses the most aggressive and $(2 .)$ uses the least aggressive stepping size, $(3 .)$ is somewhere in the middle hence the name intermediate. So from here can almost expect that if $(2 .)$ doesn’t converge then so does $(3 .)$ and $(1 .)$ . NOTE: Method $(3 .)$ always generates a monotonically decreasing sequence i.e. $f (x_{k + 1}) \leq f (x_{k})$ , and indeed this also holds for method $(2 .)$ as given in the following theorem.

Theorem 3 (Global Convergence). [1] The following holds for any $k > 0$ and let $x_{k}$ denote the newton iterates,

1.: (Standard Newton’s method) No guarantee on monotonically decreasing.
2.: (Damped Newton’s method) $f (x_{k + 1}) \leq f (x_{k}) - \frac{1}{M_{f}^{2}} ω (M_{f} λ_{f} (x_{k})),$
where $ω (x) = x - ln (1 + x)$ .
3.: (Intermediate Newton’s method) $f (x_{k + 1}) \leq f (x_{k}) - [\frac{λ_{f} {(x_{k})}^{2}}{2 \sum_{i = 0}^{2} {(M_{f} λ_{f} (x_{k}))}^{2}} + \frac{λ_{f} {(x_{k})}^{2}}{2 (1 + M_{f} λ_{f} (x_{k})) (3 + 2 M_{f} λ_{f} (x_{k})))}]$

From the above global convergence theorem we can see that if $λ_{f} (x_{k})$ is relatively large then Damped Newton’s method would be better. We have the following theorem for their convergence behavior.

Theorem 4 (Convergence Rate). [1] Let $x \in dom (f)$ and $λ = λ_{f} (x)$ , let $x^{'}$ be the point generated by one step of newton’s method (either $(1 .), (2 .), or (3 .)$ ). Then $x^{'} \in domf$ and the following holds,

1.: (Standard Newton’s method) If $λ M_{f} < 1$ then $λ_{f} (x^{'}) \leq {(\frac{M_{f} λ}{1 - M_{f} λ})}^{2} .$
2.: (Damped Newton’s method) $λ_{f} (x^{'}) \leq M_{f} λ^{2} (1 + \frac{1}{1 + M_{f} λ}) .$
3.: (Intermediate Newton’s method) If $\sum_{i = 1}^{3} {(M_{f} λ)}^{i} \leq 1$ then $λ_{f} (x^{'}) \leq M_{f} λ^{2} (1 + M_{f} λ + \frac{1 + M_{f} λ}{1 + M_{f} λ + {(M_{f} λ)}^{2}}) \leq M_{f} λ^{2} (1 + 2 M_{f} λ) .$

Notice that even though the stepping size of $(1 .)$ is more aggressive the convergence analysis doesn’t actually yield a better bound compared to $(3 .)$ or $(2 .)$ (why? consider when $λ M_{f} \approx 1$ then it yields a very bad bound), this suggests that we should completely abandon $(1 .)$ . Now method $(2 .)$ is almost a linear convergence since when $λ M_{f}$ is very big then the term $\frac{M_{f} λ^{2}}{1 + M_{f} λ}$ dominates Method $(2 .)$ does not require any condition on how close our current iterate is to the optimal solution so the idea is then always run method $(2 .)$ then switch to method $(3 .)$ when we satisfy the condition $\sum_{i = 1}^{3} {(M_{f} λ)}^{i} \leq 1$ which is yields a quadratic convergence at the end.

2 Implementation

We wish to adapt Newton’s method to any function we wish including log barrier functions and potentially many other barrier functions. Essentially Newton’s algorithm can be abstracted as follows,

Input: Initial guess

x_{0}

, function

f

, derivative

f^{'}

k \leftarrow 0

; while Not converged do

x_{k + 1} \leftarrow x_{k} + α_{k} d_{k}

;

k \leftarrow k + 1

; end Algorithm 1: Newton’s Method

Both $α_{k}, d_{k}$ may depend on the gradient and hessian. The condition to check convergence may depend on the function value, gradient, current solution, and current iterates. There are lots of ways to do this but one easy way is to: Implement a main newton routine and let the user provide lambda functions to compute the required information. But there are several caveats to this approach that I can think of,

1.: First this method does not allow compiler optimiztion in particular function inlining. To my knowledge, nor runtime polymorphism can achieve this (though I think with more aggressive compiler optimization you can still achieve this). So instead we will need compile time polymorphism (speficially CRTP Curiously Recurring Template Pattern) that uses template. Just an aside: One may ask why would function inlining be important here? Since the only difference between enabling compile time polymorphism and runtime polymorphism here is to inline the actual call to gradient/hessian and no matter what we do the calls inside of gradient/hessian would already be optimized so the only difference is that we do not need to pay the overhead of virtual table in runtime polymorphism? Not really, usually computing the gradient and hessian could involve the same intermediate computation, if we then inlined both gradient/hessian then compiler’s can avoid the same computation.
2.: If we used compile time polymorphism then usually this creates a bunch of weird error messages that are harder to resolve. To address this problem and to increase the readability of our code we use Concept which acts as a compile time contract.
3.: Recently I got pretty interested into template programming and how it can be used for compile time optimization and so I decided to implement some toy example.

1#ifndef NEWTON 
2#define NEWTON 
3 
4#include <Eigen/Dense> 
5 
6using namespace Eigen; 
7 
8template <typename T> 
9concept ObjectiveImpl = requires(const T &obj, const VectorXd &x, int it) { 
10    { obj.value_impl(x) } -> std::convertible_to<double>; 
11    { obj.gradient_impl(x) } -> std::convertible_to<VectorXd>; 
12    { obj.hessian_impl(x) } -> std::convertible_to<MatrixXd>; 
13    { obj.newtonDirection_impl(x) } -> std::same_as<VectorXd>; 
14    { obj.converged_impl(x, it) } -> std::same_as<bool>; 
15    { obj.stepsize_impl(x, it) } -> std::same_as<double>; 
16}; 
17 
18template <class Derived> class ObjectiveBase { 
19public: 
20    double value(const VectorXd& x) const { 
21        return static_cast<const Derived*>(this)->value_impl(x); 
22    } 
23 
24    VectorXd gradient(const VectorXd& x) const { 
25        return static_cast<const Derived*>(this)->gradient_impl(x); 
26    } 
27 
28    MatrixXd hessian(const VectorXd& x) const { 
29        return static_cast<const Derived*>(this)->hessian_impl(x); 
30    } 
31 
32    MatrixXd newtonDirection(const VectorXd& x) const { 
33        return static_cast<const Derived*>(this)->newtonDirection_impl(x); 
34    } 
35 
36    bool converged(const VectorXd& x, int it) const { 
37        return static_cast<const Derived*>(this)->converged_impl(x, it); 
38    } 
39 
40    double stepsize(const VectorXd& x, int it) const { 
41        return static_cast<const Derived*>(this)->stepsize_impl(x, it); 
42    } 
43}; 
44 
45template <ObjectiveImpl Obj> 
46class NewtonSolver { 
47public: 
48    NewtonSolver(Obj& objective) : obj_(objective) {} 
49 
50    bool solve() { 
51        int i = 0; 
52        while(!obj_.converged(cur_sol_, i)) { 
53            cur_sol_ = cur_sol_ + obj_.stepsize(cur_sol_, i) * obj_.newtonDirection(cur_sol_); 
54            ++i; 
55        } 
56        return true; 
57    } 
58 
59    // TODO: Currently max_iter_ and tol_ are not used... 
60    NewtonSolver& set_max_iter(int n) { max_iter_ = n; return *this; } 
61    NewtonSolver& set_tolerance_max_iter(double tol) { tol_ = tol; return *this; } 
62    NewtonSolver& set_initial_sol(const VectorXd& x) { cur_sol_ = x; return *this; } 
63    VectorXd get_sol() { return cur_sol_; } 
64private: 
65    Obj& obj_; 
66    int max_iter_ = 100; 
67    double tol_ = 1e-8; 
68    VectorXd cur_sol_; 
69}; 
70 
71#endif // NEWTON

As on how to use this library, we need to provide an derived class that inherits from ObjectiveBase that provides all the information NewtonSolver needs. We provide an example of the log barrier function,

1#ifndef LOGBARRIER 
2#define LOGBARRIER 
3 
4#include "newton.h" 
5 
6struct LogBarrier : ObjectiveBase<LogBarrier> { 
7    double value_impl(const VectorXd& x) const; 
8    VectorXd gradient_impl(const VectorXd& x) const; 
9    MatrixXd hessian_impl(const VectorXd& x) const; 
10    VectorXd newtonDirection_impl(const VectorXd& x) const; 
11    bool converged_impl(const VectorXd& x, int it) const; 
12    double stepsize_impl(const VectorXd& x, int it) const; 
13 
14    // Assume of the program min c^T x s.t. Ax <= b. 
15    LogBarrier(const MatrixXd& A, const VectorXd& b, const VectorXd& c, double t = 1) : A_(A), b_(b), c_(c), t_(t), m_(A.rows()), n_(A.cols()) {} 
16 
17private: 
18    MatrixXd A_; 
19    VectorXd b_; 
20    VectorXd c_; 
21    double t_; 
22    int m_, n_; 
23}; 
24 
25#endif // LOGBARRIER

1#include "logbarrier.h" 
2 
3double LogBarrier::value_impl(const VectorXd& x) const { 
4    return t_ * c_.dot(x) - (b_ - A_ * x).array().log().sum(); 
5} 
6 
7VectorXd LogBarrier::gradient_impl(const VectorXd& x) const { 
8    VectorXd s = b_ - A_ * x; 
9    return t_ * c_ + A_.transpose() * s.cwiseInverse(); 
10} 
11 
12MatrixXd LogBarrier::hessian_impl(const VectorXd& x) const { 
13    VectorXd s = b_ - A_ * x; 
14    VectorXd s_inv_sq = s.array().square().inverse(); 
15    return A_.transpose() * s_inv_sq.asDiagonal() * A_; 
16} 
17 
18VectorXd LogBarrier::newtonDirection_impl(const VectorXd& x) const { 
19    return -(hessian_impl(x).lu().solve(gradient_impl(x))); 
20} 
21 
22bool LogBarrier::converged_impl(const VectorXd& x, int it) const { 
23    return it >= 27; 
24} 
25 
26double LogBarrier::stepsize_impl(const VectorXd& x, int it) const { 
27    // TODO: Should optimize this to resuse previously computed values if possible. 
28    double lambda = sqrt(gradient_impl(x).dot(hessian_impl(x).lu().solve(gradient_impl(x)))); 
29    double alpha; 
30 
31    double m_lambda = m_*lambda; 
32 
33    if (m_lambda <= 0.5) { 
34        alpha = (1+m_lambda)/(1+m_lambda+m_lambda*m_lambda); 
35    } else { 
36        alpha = 1/(1+m_lambda); 
37    } 
38 
39    return alpha; 
40}

Further detailed code can be found on my github page.

References

[1] Yurii Nesterov. Lectures on Convex Optimization, volume 137 of Springer Optimization and Its Applications. Springer, 2 edition, 2018.