Robotics - Gauss Newton (GN) and Levenberg Marquardt (LM) Optimizers - Rico贾若童的博客

Newton’s Method

Solving An Equation

To find an arbitrary equation’s root $f (x) = 0$ ,

We start from an arbitrary point $x_{0}$ that’s hopefully close to the solution, $x_{s}$
The main idea of Newton’s method is, we draw a tangent line at $x_{0}$ , this line will cross y=0 at $x_{1}$ . It’s most likely that this point is closer to $x_{s}$ than $x_{0}$ . So, this is to solve $x_{1} f (x_{0}) - x_{0} f^{'} (x_{0}) + f (x_{0}) = 0$ , and we get $x_{1} = x_{0} - \frac{f (x_{0})}{f^{'} (x_{0})}$

So at each iteration, the newer estimate is:

\begin{matrix} \begin{aligned} x_{n + 1} = x_{n} - \frac{f (x_{0})}{f^{'} (x_{0})} \end{aligned} \end{matrix}

Some caution-worthy notes are:

$f^{'} (x_{0})$ should be non-zero. Otherwise, the estimate will NOT move
In vanilla gradient descent, we only find the gradient and issue an update $Δ x = λ g$ using a fixed step size $λ$ . However, in Newton’s method, we not only follow the gradient direction, but we also “solve” for the step size. That makes Newton’s method faster.

Example: Solve For Square-Roots

If we want find $\sqrt{m}$ , we can format this problem to solving:

\begin{matrix} \begin{aligned} f (x) = x^{2} - m = 0 \\ \to x_{n + 1} = x_{n} - \frac{x_{n}^{2} - m}{2 x_{n}} = \frac{x_{n}}{2} + \frac{m}{2 x_{n}} \end{aligned} \end{matrix}

const double EPS = 0.0001; 
double sqrt(double m){
    if (m == 0) return m;
    double x = m;
    double last_res;
    do {
        last_res = x;
        x = x / 2.0 + m / (2 * x);
    } while (abs(last_res - x) > EPS);
    return x;
}

Faster version (TODO), inspired by John Carmack

double sqrt_c(float x) { 
    if(x == 0) return 0; 
    float result = x; 
    float xhalf = 0.5f*result; 
    int i = *(int*)&result; 
    // 1048009350
    i = 0x5f375a86- (i>>1); // what the fuck? 
    result = *(float*)&i; 
        cout<<result<<endl; //0.241556
    result = result*(1.5f-xhalf*result*result); // Newton step, repeating increases accuracy 
    result = result*(1.5f-xhalf*result*result); 
    return 1.0f/result; 
}

Newton’s Method for Optimization

When optimizing f(x), again we start from an arbitrary point $x_{0}$ . If we can achieve the optimum at $x_{0} + Δ x$ :

We think of it as:

\begin{matrix} \begin{aligned} Taylor Expansion: \\ f (x + Δ x) \approx f (x_{0}) + f^{'} (x_{0}) Δ x + \frac{1}{2} (f^{″} (x_{0}))^{2} Δ x^{2} \\ Optimum: \\ \frac{\partial f (x_{0} + Δ x)}{\partial Δ x} = 0 \\ \to f (x + Δ x)^{'} = f^{'} (x_{0}) + f^{″} (x_{0}) Δ x = 0 \\ \to Δ x = - \frac{f^{'} (x_{0})}{f^{″} (x_{0})} \\ \to x_{1} = x_{0} + Δ x = x_{0} - \frac{f^{'} (x_{0})}{f^{″} (x_{0})} \end{aligned} \end{matrix}

So iteratively, at timestep n, we have:

\begin{matrix} \begin{aligned} x_{n + 1} = x_{n} - \frac{f^{'} (x_{n})}{f^{″} (x_{n})} \end{aligned} \end{matrix}

See how Newton’s method has similar forms for both equation solving and optimization? The better the function f is approximated by Taylor Expansion, the faster we converge.

For higher dimensions, if f is a scalar-valued function, and x is [m,1] we have:

\begin{matrix} \begin{aligned} x_{n + 1} = x_{n} - J (x_{n}) H (x_{n})^{- 1} \end{aligned} \end{matrix}

J is [1, m], H is [m,m]

Gauss-Newton Optimization

In Gauss Newton, we specifically look at minimizing a least squares problem. Assume we have a:

scalar-valued cost function $c (x)$ ,
vector-valued function: $f (x)$ , [m, 1]
Jacobian $J_{0}$ at $x_{0}$ is consequently [m, n]
Hessian $H$ is $D^{2} c (x)$ . It’s approximated as $J^{T} J$

\begin{matrix} \begin{aligned} c (x) = | f (x)^{2} | \\ x * = a r g m i n (| f (x)^{2} |) \\ First order Taylor Expansion: \\ a r g m i n_{Δ x} (| f (x + Δ x)^{2} |) \\ = a r g m i n_{Δ x} [(f (x_{0}) + J_{0} Δ x)^{T} (f (x_{0}) + J_{0} Δ x)] \\ = a r g m i n_{Δ x} [f (x_{0})^{T} f (x_{0}) + f (x_{0})^{T} J_{0} Δ x + (J_{0} Δ x)^{T} f (x_{0}) + (J_{0} Δ x)^{T} (J_{0} Δ x)] \\ = a r g m i n_{Δ x} [f (x_{0})^{T} f (x_{0}) + 2 f (x_{0})^{T} J_{0} Δ x + (J_{0} Δ x)^{T} (J_{0} Δ x)] \end{aligned} \end{matrix}

Take the derivative of the above and set it to 0, we get

\begin{matrix} \begin{aligned} \frac{\partial f (x + Δ x)^{2}}{\partial Δ x} = 2 J_{0}^{T} f (x_{0}) + [(J_{0}^{T} J_{0}) + (J_{0}^{T} J_{0})^{T}] Δ x \\ = 2 J_{0}^{T} f (x_{0}) + 2 (J_{0}^{T} J_{0}) Δ x \\ = 0 \end{aligned} \end{matrix}

So we can solve for $Δ x$ with $H = J_{0}^{T} J_{0}$ , $b = - J_{0}^{T} f (x_{0})$ :

\begin{matrix} \begin{aligned} (J_{0}^{T} J_{0}) Δ x = - J_{0}^{T} f (x_{0}) \\ \to H Δ x = g \end{aligned} \end{matrix}

Note: because $J_{0}$ may not have an inverse, here we cannot multiply $J_{0}^{- 1}$ to eliminate $J_{0}^{T}$
In fact, to $Δ x$ is available if and only if $H$ is positive definite.
In least square, $f (x)$ is a.k.a residuals. Usually, it represents the error between a data point and from its ground truth.

In SLAM, we always frame this least squares problem with e = [observered_landmark - predicted_landmark] at each landmark. So all together, we want to minimize the total least squares of the difference between observations and predictions. In the meantime, at each landmark, there is an error covariance, so all together, there’s an error matrix $Σ$ . Here in cost calculation, we take $Σ^{- 1}$ so the larger the error covariance, the lower the weight the corresponding difference gets.

With $e (x + Δ x) \approx e (x) + J Δ x$ ,

\begin{matrix} \begin{aligned} x * = a r g m i n (| e^{T} Σ^{- 1} e |) \\ \to a r g m i n_{Δ x} (| e (x + Δ x)^{T} Σ^{- 1} e (x + Δ x) |) \\ similar steps as above ... \\ \to (J_{0}^{T} Σ^{- 1} J_{0}) Δ x = - J_{0}^{T} Σ^{- 1} f (x_{0}) \end{aligned} \end{matrix}

Using Cholesky Decomposition, one can get $Σ^{- 1} = A^{T} A$ . Then we can write the above as

\begin{matrix} \begin{aligned} ((A J_{0})^{T} (A J_{0})) Δ x = - (A J_{0})^{T} A f (x_{0}) \end{aligned} \end{matrix}

For a more detailed derivation, please see here

When $Δ x$ is small, or $f (x)$ is small, GN has faster convergence. When either or both of them are large, or when the target function is highly non-linear, G-N does not work well.

Levenberg-Marquardt (LM) Optimization

Again, Taylor expansion works better when $Δ x$ is small, so the function can be better estimated by it. So, similar to regularization techniques on step sizes in deep learning, like L1, L2 regularization, we can regularize the step size, $Δ x$ . $μ$ is called “the damping factor”

\begin{matrix} \begin{aligned} (H + μ I_{H}) Δ x = - J_{0}^{T} f (x_{0}) = g \end{aligned} \end{matrix}

Intuitively,

as $μ$ grows, the diagonal identity matrix $μ I_{H}$ grows, so $H + μ I_{H} \to μ I_{H}$ . So, $Δ x \approx (H + μ I_{H})^{- 1} g = \frac{g}{μ}$ , which means $Δ x$ grows smaller. In the meantime, $Δ x$ will be similar to that in gradient descent.
as $μ$ becomes smaller, $Δ x$ will become more like Gauss-Newton. However, due to $μ I_{H}$ , $(H + μ I_{H})$ is positive semi-definite, which provides more stability for solving for $Δ x$ .

LM in general is more stable than GN. However, it’s more computationally costly.

References

https://scm_mos.gitlab.io/algorithm/newton-and-gauss-newton/

Related Issues not found

Please contact @ricojia to initialize the comment

Robotics - Gauss Newton (GN) and Levenberg Marquardt (LM) Optimizers

Newton's Method, GN, LM optimizers