8. Quasi-Newton-Raphson Methods by MIT OpenCourseWare

Summary by www.lecturesummary.com: 8. Quasi-Newton-Raphson Methods by MIT OpenCourseWare

• Shift to Nonlinear Equations and Iterative Techniques (0:00 - 1:06)

  • The lecture shifts from solving systems of linear equations to solving systems of nonlinear equations.
  • Linear algebra is at the center of methods applied to nonlinear equations.
  • Iterative methods are required since these problems are complex; the number and location of solutions are unclear, and no precise methods of locating them are available.
  • Iterative methods reduce nonlinear equations into tractable problems, namely iterates of systems of linear equations.
  • The primary method outlined is the Newton-Raphson method.
  • The lecture will discuss techniques, known as quasi-Newton-Raphson techniques, for overcoming challenges and limitations of the basic Newton-Raphson technique.
  • These techniques are selected as required when the Newton-Raphson technique poses challenges.
  • Programs such as MATLAB's fsolve utilize elements of these quasi-Newton-Raphson techniques.

  Scale of Nonlinear Systems (1:06 - 2:08)

  • It is possible to have systems of nonlinear equations of arbitrarily large sizes.
  • A good example is the steady Navier-Stokes equations, which is a nonlinear partial differential equation of fluid velocity and pressure.
  • The discretization of such equations results in a system of coupled nonlinear equations at every point in the fluid of interest.
  • The amount of points, and hence equations, may be enormous.
  • Even though fluid mechanics has other approaches at times, theoretically, approaches such as these must be used in order to solve systems of any quantity of nonlinear equations.

  Newton-Raphson Method Summary and Babylonian Method Explanation (2:08 - 4:53)

  • The Newton-Raphson technique relies on the principle of linearization.
  • Starting from a preliminary guess for the solution, the technique linearizes the function around the guess and determines where the linearized version has an intercept.
  • This intercept is taken as the next, and hopefully better, guess, and the process is repeated.
  • The Newton-Raphson method is locally convergent, i.e., if the initial estimate is close enough to the root, the process is guaranteed to converge to the root.
  • The Babylonian method for obtaining square roots can be deduced from the 1D Newton-Raphson formula.
  • The derivative of f(x) = x² - S is f'(x) = 2x.
  • Replacing f(x) and f'(x) in the 1D Newton-Raphson formula (xᵢ₊₁ = xᵢ - f(xᵢ)/f'(xᵢ)) gives the Babylonian method formula.
  • The Babylonian method has quadratic convergence to the square root.
  • Good initial guesses are crucial since distances far from the solution need more iterations to get there.

  Graphical View of Newton-Raphson in Multi-Dimensions (4:53 - 6:17)

  • In several dimensions, for a function with several components (e.g., f₁, f₂) and several unknowns (e.g., x₁, x₂), the derivative is substituted by the Jacobian matrix.
  • Linearization of the function at a guess is not a line but a plane (or hyperplane).
  • This plane crosses the x₁x₂ plane (or n-dimensional unknown space) in a line.
  • For a system of equations, linearizations of each component function give lines (or hyperplanes) in the unknown space.
  • The intersection of the lines (or hyperplanes) provides the next best estimate of the solution.
  • To find this intersection, calculate Jacobian inverse times f.
  • Mapping to the space of unknowns, the iterates trace a path that converges to the locally unique solution ultimately.
  • Beginning near the root results in rapid convergence.

  • Quadratic Convergence Proof (6:17 - 7:56)

    The Newton-Raphson method is quadratically convergent.

    • In 1D, the step i+1 error (deviation from true root x*) is proportional to the step i error by |xᵢ₊₁ - x*| ≈ ½ |f''(x*)/f'(x*)| * |xᵢ - x*|².
    • As i approaches infinity, the ratio of the absolute error in step i+1 to the absolute error in step i squared is bounded above by a constant.
    • This condition |errorᵢ₊₁| / |errorᵢ|² ≤ Constant is the definition of quadratic convergence.
    • Quadratic convergence implies that the number of correct digits roughly doubles on every iteration.
    • Provided that the derivative f'(x*) at the root is not zero, this is the case.
    • If f'(x*) is zero, then the analysis is invalid, and quadratic convergence is lost.
    • In the multi-dimensional situation, the determinant of the Jacobian at the root takes the place of the derivative.
    • Quadratic convergence is lost if the Jacobian is singular at the root, leading to linear convergence instead.
    • Quadratic convergence is only assured if the initial guess is sufficiently close to the root. With poor initial guesses, convergence is not assured, and the method could diverge.

    Problems with the Newton-Raphson Method (7:56 - 9:50)

    • Hung up: The algorithm may oscillate between local minima or maxima, with the linearization pointing away from the root in the next step. It is locally convergent but not globally.
    • Divergence: Iterates can blow up along asymptotes.
    • Overshooting: Iterations can continuously overshoot the root, often resulting in divergence or slow convergence at roots where the derivative does not exist.
    • Basins of Attraction: In the case of functions with multiple roots, convergence to different roots is achieved by using different initial guesses. The areas in the unknown space from where initial guesses converge to a specific root are referred to as basins of attraction.
    • Such basins may be fractal in character, making it hard to anticipate which root the method will converge to from a given initial guess. Even a slight variation of the initial guess could lead to convergence to another root.
    • One often doesn't know how many roots there are or which one will be found; the found root might be unphysical or not the desired one. This is a problem, even for relatively simple functions. Polynomial equations are particularly prone to this.

    Further Weaknesses and Need for Modifications (9:50 - 10:59)

    • Difficulty of Analytic Jacobian Calculation: It can be cumbersome, error-ridden, or even impossible to compute the Jacobian analytically for complex functions or those including simulations or data.
    • Complexity of Jacobian Inversion: It is computationally costly to invert the Jacobian matrix, particularly for large-scale systems. This linear system solving at every iteration is undesirable.
    • Failure of Global Convergence/Convergence to Nearest Root: The method could fail to converge at all or fail to converge to the nearest root because of overshoot or the basin of attraction.
    • Quasi-Newton-Raphson algorithms adjust the linearization to help mitigate these problems.
    • These modifications come with a penalty: the loss of quadratic convergence, resulting in slower convergence rates. However, they can improve reliability or make individual iterations computationally cheaper, potentially leading to faster overall calculation.

    Approximating the Jacobian with Finite Differences (10:59 - 13:15)

    • One modification is to approximate the Jacobian numerically using finite differences.
    • This is required when there is no analytical formula for the function or its derivatives available, e.g., when function values are available from simulations or data interpolation.
    • The partial derivative ∂fᵢ/∂xⱼ can be approximated by (fᵢ(x + εeⱼ) - fᵢ(x)) / ε, where eⱼ is a unit vector in the j-th direction.
    • Selecting the finite difference step size ε is critical; a rule of thumb is to take ε ≈ √(machine precision) * |x|. Very small ε can cause larger error because of floating-point truncation errors.
    • Calculating the entire Jacobian matrix at a point by finite differences takes n+1 function evaluations for an n-equation system. This is still computationally costly if function evaluations are expensive (e.g., from simulations).
    • Approximating the Jacobian lowers the rate of convergence from quadratic to superlinear. The rate is a function of how good the approximation is and how sensitive the function is.

    • The fsolve function in MATLAB applies finite difference approximations of the Jacobian as a fallback if an analytical Jacobian is not provided.

      Secant Method and Broden's Method (13:15 - 16:21)

      The secant method in 1D employs a less accurate estimate for the derivative:

      • The derivative at xᵢ is estimated as (f(xᵢ) - f(xᵢ₋₁)) / (xᵢ - xᵢ₋₁).
      • This method cannot be straightforwardly extended to more than one dimension.
      • The multi-dimensional counterpart of the secant equation is J * (xᵢ - xᵢ₋₁) = f(xᵢ) - f(xᵢ₋₁).
      • This equation is underdetermined for the N² elements of the Jacobian with just the N-dimensional step vector and function difference vector.
      • Broden's approach is a method for selecting a good approximation to the Jacobian from among the numerous solutions to this underdetermined problem.
      • It updates the Jacobian approximation iteratively by a rank-one update formula.
      • One of Broden's advantages is that the inverse of the new Jacobian Jᵢ can be obtained directly from the inverse of the old Jacobian Jᵢ₋₁ by applying the Sherman-Morrison formula.
      • This skips the costly operation of solving a linear system of equations at every iteration; rather, the Jacobian inverse matrix is iteratively updated.
      • Precision is sacrificed, but computation time for each iteration is saved.

      Damped Newton-Raphson Method (16:21 - 22:56)

      The Newton-Raphson step in standard form sometimes lands on a new point where the function value is larger than at the current point:

      • This is typical far from the root.
      • The Newton-Raphson step gives a direction reducing the function value locally, but the overall step size may overshoot.
      • The damped Newton-Raphson method adds a damping factor α (between 0 and 1) to adjust the step: xᵢ₊₁ = xᵢ - α * J⁻¹(xᵢ) * f(xᵢ).
      • The intention is to select α so that the function norm at the new point ||f(xᵢ₊₁)|| is minimized.
      • There are practical approaches based on approximate methods to select α.
      • One popular method is a line search, e.g., the Armijo line search.
      • Armijo line search begins with α=1 (the entire Newton-Raphson step).
      • It then tests whether the norm of the function at the new point ||f(xᵢ₊₁)|| is suitably smaller than the norm at the current point ||f(xᵢ)||.
      • If the step is good (e.g., ||f(xᵢ₊₁)|| < ||f(xᵢ)||), α=1 is accepted, and the step is executed.
      • Otherwise, α is decreased (e.g., set to α/2), and the test is tried again with the shorter step.
      • This is repeated until a step size is identified which decreases the function norm, and this step is then taken.
      • Since the Newton-Raphson direction ensures local reduction in function value if steps are small enough, there will always be a suitable α (unless a minimum/maximum has been reached).
      • The Newton-Raphson method with damping is globally convergent, in the sense that it is guaranteed to converge eventually.
      • But it is not globally convergent to a root; instead, it will converge to local minima or maxima of the function norm.
      • Damping reduces the basins of attraction, reducing the convergence behavior to potentially being less fractal and more deterministic.
      • Near a root, α=1 is generally accepted straight away, and the technique regains its fast convergence.
      • The additional expense of the damping/line search phase per iteration is comparatively inexpensive relative to the computation or inversion of the Jacobian.
      • MATLAB's fsolve tends to employ damping or step sizing limitation automatically.
      • One can use a matrix instead of the scalar α for damping to alter the direction or scale differently in varying dimensions.
      • The cost of damping is a reduced rate of convergence but may enhance robustness and reliability.