parameters: (n, m, uf, s, V, Delta, initial_alpha=None,
                             rtol=0.01, max_iter=10)
return_statement: return p, alpha, it + 1
docstring: Solve a trust-region problem arising in least-squares minimization.
      
      This function implements a method described by J. J. More [1]_ and used
      in MINPACK, but it relies on a single SVD of Jacobian instead of series
      of Cholesky decompositions. Before running this function, compute:
      ``U, s, VT = svd(J, full_matrices=False)``.
      
      Parameters
      ----------
      n : int
          Number of variables.
      m : int
          Number of residuals.
      uf : ndarray
          Computed as U.T.dot(f).
      s : ndarray
          Singular values of J.
      V : ndarray
          Transpose of VT.
      Delta : float
          Radius of a trust region.
      initial_alpha : float, optional
          Initial guess for alpha, which might be available from a previous
          iteration. If None, determined automatically.
      rtol : float, optional
          Stopping tolerance for the root-finding procedure. Namely, the
          solution ``p`` will satisfy ``abs(norm(p) - Delta) < rtol * Delta``.
      max_iter : int, optional
          Maximum allowed number of iterations for the root-finding procedure.
      
      Returns
      -------
      p : ndarray, shape (n,)
          Found solution of a trust-region problem.
      alpha : float
          Positive value such that (J.T*J + alpha*I)*p = -J.T*f.
          Sometimes called Levenberg-Marquardt parameter.
      n_iter : int
          Number of iterations made by root-finding procedure. Zero means
          that Gauss-Newton step was selected as the solution.
      
      References
      ----------
      .. [1] More, J. J., "The Levenberg-Marquardt Algorithm: Implementation
             and Theory," Numerical Analysis, ed. G. A. Watson, Lecture Notes
             in Mathematics 630, Springer Verlag, pp. 105-116, 1977.
docstring_summary: Solve a trust-region problem arising in least-squares minimization.
      
      This function implements a method described by J. J. More [1]_ and used
      in MINPACK, but it relies on a single SVD of Jacobian instead of series
      of Cholesky decompositions. Before running this function, compute:
      ``U, s, VT = svd(J, full_matrices=False)``.
      
      Parameters
      ----------
      n : int
          Number of variables.
      m : int
          Number of residuals.
      uf : ndarray
          Computed as U.T.dot(f).
      s : ndarray
          Singular values of J.
      V : ndarray
          Transpose of VT.
      Delta : float
          Radius of a trust region.
      initial_alpha : float, optional
          Initial guess for alpha, which might be available from a previous
          iteration. If None, determined automatically.
      rtol : float, optional
          Stopping tolerance for the root-finding procedure. Namely, the
          solution ``p`` will satisfy ``abs(norm(p) - Delta) < rtol * Delta``.
      max_iter : int, optional
          Maximum allowed number of iterations for the root-finding procedure.
      
      Returns
      -------
      p : ndarray, shape (n,)
          Found solution of a trust-region problem.
      alpha : float
          Positive value such that (J.T*J + alpha*I)*p = -J.T*f.
          Sometimes called Levenberg-Marquardt parameter.
      n_iter : int
          Number of iterations made by root-finding procedure. Zero means
          that Gauss-Newton step was selected as the solution.
      
      References
      ----------
      .. [1] More, J. J., "The Levenberg-Marquardt Algorithm: Implementation
             and Theory," Numerical Analysis, ed. G. A. Watson, Lecture Notes
             in Mathematics 630, Springer Verlag, pp. 105-116, 1977.
docstring_tokens: Solve a trust - region problem arising in least - squares minimization . This function implements a method described by J . J . More [ 1 ] _ and used in MINPACK but it relies on a single SVD of Jacobian instead of series of Cholesky decompositions . Before running this function compute : U s VT = svd ( J full_matrices = False ) . Parameters ---------- n : int Number of variables . m : int Number of residuals . uf : ndarray Computed as U . T . dot ( f ) . s : ndarray Singular values of J . V : ndarray Transpose of VT . Delta : float Radius of a trust region . initial_alpha : float optional Initial guess for alpha which might be available from a previous iteration . If None determined automatically . rtol : float optional Stopping tolerance for the root - finding procedure . Namely the solution p will satisfy abs ( norm ( p ) - Delta ) < rtol * Delta . max_iter : int optional Maximum allowed number of iterations for the root - finding procedure . Returns ------- p : ndarray shape ( n ) Found solution of a trust - region problem . alpha : float Positive value such that ( J . T * J + alpha * I ) * p = - J . T * f . Sometimes called Levenberg - Marquardt parameter . n_iter : int Number of iterations made by root - finding procedure . Zero means that Gauss - Newton step was selected as the solution . References ---------- .. [ 1 ] More J . J . The Levenberg - Marquardt Algorithm : Implementation and Theory Numerical Analysis ed . G . A . Watson Lecture Notes in Mathematics 630 Springer Verlag pp . 105 - 116 1977 .
function: def solve_lsq_trust_region(n, m, uf, s, V, Delta, initial_alpha=None,
                             rtol=0.01, max_iter=10):
      """Solve a trust-region problem arising in least-squares minimization.
      
      This function implements a method described by J. J. More [1]_ and used
      in MINPACK, but it relies on a single SVD of Jacobian instead of series
      of Cholesky decompositions. Before running this function, compute:
      ``U, s, VT = svd(J, full_matrices=False)``.
      
      Parameters
      ----------
      n : int
          Number of variables.
      m : int
          Number of residuals.
      uf : ndarray
          Computed as U.T.dot(f).
      s : ndarray
          Singular values of J.
      V : ndarray
          Transpose of VT.
      Delta : float
          Radius of a trust region.
      initial_alpha : float, optional
          Initial guess for alpha, which might be available from a previous
          iteration. If None, determined automatically.
      rtol : float, optional
          Stopping tolerance for the root-finding procedure. Namely, the
          solution ``p`` will satisfy ``abs(norm(p) - Delta) < rtol * Delta``.
      max_iter : int, optional
          Maximum allowed number of iterations for the root-finding procedure.
      
      Returns
      -------
      p : ndarray, shape (n,)
          Found solution of a trust-region problem.
      alpha : float
          Positive value such that (J.T*J + alpha*I)*p = -J.T*f.
          Sometimes called Levenberg-Marquardt parameter.
      n_iter : int
          Number of iterations made by root-finding procedure. Zero means
          that Gauss-Newton step was selected as the solution.
      
      References
      ----------
      .. [1] More, J. J., "The Levenberg-Marquardt Algorithm: Implementation
             and Theory," Numerical Analysis, ed. G. A. Watson, Lecture Notes
             in Mathematics 630, Springer Verlag, pp. 105-116, 1977.
      """
      def phi_and_derivative(alpha, suf, s, Delta):
          """Function of which to find zero.
          
          It is defined as "norm of regularized (by alpha) least-squares
          solution minus `Delta`". Refer to [1]_.
          """
          denom = s**2 + alpha
          p_norm = norm(suf / denom)
          phi = p_norm - Delta
          phi_prime = -np.sum(suf ** 2 / denom**3) / p_norm
          return phi, phi_prime
  
      suf = s * uf
  
      # Check if J has full rank and try Gauss-Newton step.
      if m >= n:
          threshold = EPS * m * s[0]
          full_rank = s[-1] > threshold
      else:
          full_rank = False
  
      if full_rank:
          p = -V.dot(uf / s)
          if norm(p) <= Delta:
              return p, 0.0, 0
  
      alpha_upper = norm(suf) / Delta
  
      if full_rank:
          phi, phi_prime = phi_and_derivative(0.0, suf, s, Delta)
          alpha_lower = -phi / phi_prime
      else:
          alpha_lower = 0.0
  
      if initial_alpha is None or not full_rank and initial_alpha == 0:
          alpha = max(0.001 * alpha_upper, (alpha_lower * alpha_upper)**0.5)
      else:
          alpha = initial_alpha
  
      for it in range(max_iter):
          if alpha < alpha_lower or alpha > alpha_upper:
              alpha = max(0.001 * alpha_upper, (alpha_lower * alpha_upper)**0.5)
  
          phi, phi_prime = phi_and_derivative(alpha, suf, s, Delta)
  
          if phi < 0:
              alpha_upper = alpha
  
          ratio = phi / phi_prime
          alpha_lower = max(alpha_lower, alpha - ratio)
          alpha -= (phi + Delta) * ratio / Delta
  
          if np.abs(phi) < rtol * Delta:
              break
  
      p = -V.dot(suf / (s**2 + alpha))
  
      # Make the norm of p equal to Delta, p is changed only slightly during
      # this. It is done to prevent p lie outside the trust region (which can
      # cause problems later).
      p *= Delta / norm(p)
  
      return p, alpha, it + 1
function_tokens: def solve_lsq_trust_region ( n , m , uf , s , V , Delta , initial_alpha = None , rtol = 0.01 , max_iter = 10 ) : def phi_and_derivative ( alpha , suf , s , Delta ) : """Function of which to find zero.
          
          It is defined as "norm of regularized (by alpha) least-squares
          solution minus `Delta`". Refer to [1]_.
          """ denom = s ** 2 + alpha p_norm = norm ( suf / denom ) phi = p_norm - Delta phi_prime = - np . sum ( suf ** 2 / denom ** 3 ) / p_norm return phi , phi_prime suf = s * uf # Check if J has full rank and try Gauss-Newton step. if m >= n : threshold = EPS * m * s [ 0 ] full_rank = s [ - 1 ] > threshold else : full_rank = False if full_rank : p = - V . dot ( uf / s ) if norm ( p ) <= Delta : return p , 0.0 , 0 alpha_upper = norm ( suf ) / Delta if full_rank : phi , phi_prime = phi_and_derivative ( 0.0 , suf , s , Delta ) alpha_lower = - phi / phi_prime else : alpha_lower = 0.0 if initial_alpha is None or not full_rank and initial_alpha == 0 : alpha = max ( 0.001 * alpha_upper , ( alpha_lower * alpha_upper ) ** 0.5 ) else : alpha = initial_alpha for it in range ( max_iter ) : if alpha < alpha_lower or alpha > alpha_upper : alpha = max ( 0.001 * alpha_upper , ( alpha_lower * alpha_upper ) ** 0.5 ) phi , phi_prime = phi_and_derivative ( alpha , suf , s , Delta ) if phi < 0 : alpha_upper = alpha ratio = phi / phi_prime alpha_lower = max ( alpha_lower , alpha - ratio ) alpha -= ( phi + Delta ) * ratio / Delta if np . abs ( phi ) < rtol * Delta : break p = - V . dot ( suf / ( s ** 2 + alpha ) ) # Make the norm of p equal to Delta, p is changed only slightly during # this. It is done to prevent p lie outside the trust region (which can # cause problems later). p *= Delta / norm ( p ) return p , alpha , it + 1