NonlinearLeastSquares.py
 
Version: 2.0.4
   
Author: Avinash Kak (kak@purdue.edu)
 
Date: 2026-October-12
 
 

Download Version 2.0.4:  gztar	Total number of downloads (all versions):      This count is automatically updated at every rotation of      the weblogs (normally once every two to four days)      Last updated:

View the main module file in your browser  
View the code for OptimizedSurfaceFit for optimized surface fitting to noisy height data
View the code for ProjectiveCamera for camera modeling
 
 
 
CHANGES:

  Version 2.0.4:
 
    The primary reason for this upgrade was that the foundational Scipy library no
    longer supports the trigonometric functions needed for camera modeling.  Those
    must now be accessed directly through the Numpy library.  While I was making
    those changes to the ProjectiveCamera class in the NonlinearLeastSquares module I
    discovered a couple of bugs that have now been fixed.
 
  Version 2.0.3:
 
    The previous version had two bugs in the iterative logic of the
    Levenberg-Marquardt algorithm, one in updating the current prediction error and
    the other in updating the damping coefficient.  These have been fixed in this
    version.
 
  Version 2.0.2:
 
    This version includes additional functionality in the ProjectiveCamera class
    for constructing flyby camera-image sequences of virtual objects in the world
    coordinate frame.  See the script flyby_photos_of_3D_planar_shape.py in the
    ExamplesStructureFromCameraMotion directory of the distribution for how you
    can generate flyby images of a planar object in the YZ-plane.
 
  Version 2.0.1:
 
    Cleaned up the code to make it work with the current version of the matplotlib
    library that is used for displaying the results.  I have also removed the use
    of the "is" and "is not" comparison operators for numeric and string literals
    since such use has now been deprecated.  Also added additional comment blocks
    to the main functions in the module file.
 
  Version 2.0.0:
 
    This version includes sparse bundle adjustment (SBA) to speed up code
    execution when using nonlinear least-squares for solving
    structure-from-camera-motion problems with uncalibrated cameras.  SBA exploits
    the sparseness of the Jacobian encountered in such problems.  This version
    also includes improvements to the Levenberg-Marquardt code itself. The most
    important improvement consists of how the "mu" factor (renamed "alambda" in
    this module) is initialized.  Version 2.0.0 also provides additional methods
    for the visualization of the results.
 
  Version 1.5.0:
 
    This version adds a new class, ProjectiveCamera, to the module for
    demonstrating how nonlinear least-squares can be used for estimating structure
    of a scene from the data recorded with a calibrated camera in motion.  For a
    simulated demonstration, you can create calibrated cameras from a generic
    instance of the ProjectiveCamera class that is then subject to various
    translational and rotational transformations.
 
  Version 1.1.1:
 
    This version includes constructor options to control the size and the position of
    the graphics for displaying the results of nonlinear least-squares optimization.
 
  Version 1.1:
 
    This version fixes a bug in the synthetic data generator function used for
    illustrating the functionality of the NonlinearLeastSquares class.  Changes also
    made to the information that is printed out when the module is run in the debug
    mode.
 
  Version 1.0:
 
    In the form of a class named NonlinearLeastSquares, this module provides a domain
    agnostic implementation of nonlinear least-squares algorithms (gradient-descent
    and Levenberg-Marquardt) for fitting a model to observed data.  Typically, the
    model involves several parameters and each observed data element can be expressed
    as a function of those parameters plus noise.  The goal of nonlinear
    least-squares is to estimate the best values for the parameters given all of the
    observed data.  In order to illustrate how to use the NonlinearLeastSquares
    class, the module also comes with another class, OptimizedSurfaceFit, whose job
    is to fit the best surface (of a specified analytical form) to noisy height data
    over the XY-plane. The model in this case is the analytical description of the
    surface and the goal of nonlinear least-squares is to estimate the best values
    for the parameters in the analytical description.
 
 
USAGE FOR OPTIMAL SURFACE FITTING:
 
    The module includes a domain specific class, OptimizedSurfaceFit, for
    demonstrating how nonlinear least-squares in the form of the Levenberg-Marquardt
    algorithm can be used for optimally fitting analytically specified surfaces to
    noisy data over a plane.
    
    For surface fitting, you need to create an instance of the NonlinearLeastSquares
    class as shown in line (A) below.  In addition, you need to create an instance of
    a domain-specific class such as OptimizedSurfaceFit as is shown at line (B)
    below.  Note that in the example shown in line (B), the domain-specific class has
    two duties: (1) to synthetically generate noisy height data; and (2) to call on
    NonlinearLeastSquares to fit an optimized model to the synthetically generated
    data.
 
        optimizer =  NonlinearLeastSquares(                              #(A)
                         max_iterations = 400,
                         delta_for_jacobian = 0.000001,
                         delta_for_step_size = 0.0001,        # only needed for GD
                     )
 
        surface_fitter = OptimizedSurfaceFit(                            #(B)                
                                gen_data_synthetically = True,
                                datagen_functional = "7.8*(x - 0.5)**3 + 2.2*(y - 0.5)**2",
                                data_array_size = (16,16),
                                how_much_noise_for_synthetic_data = 0.7,
                                model_functional = "a*(x-b)**3 + c*(y-d)**2",
                                initial_param_values = {'a':2.0, 'b':0.4, 'c':0.8, 'd':0.4},
                                display_needed = True,
                         )
 
        surface_fitter.set_constructor_options_for_optimizer(optimizer)  #(C)
 
        result = surface_fitter.calculate_best_fitting_surface('lm')     #(D)
 
        or 
 
        result = surface_fitter.calculate_best_fitting_surface('gd')     #(E)
 
    In line (C), it is the job of the OptimizedSurfaceFit instance as constructed in
    line (B) to export the domain-related information (after it is packaged under the
    hood into a domain-agnostic form) to the instance of NonlinearLeastSquares that
    was constructed in line (A).  The information that the OptimizedSurfaceFit
    instance conveys to the NonlinearLeastSquares instance includes:
 
        1) The observed noisy data vector.  This vector is denoted X in the
           NonlinearLeastSquares class.  The same notation is used for the
           observed data vector in the usage-example class OptimizedSurfaceFit.
 
        2) The initial values to be used for the parameters of the model.
 
        3) A vector of the predicted values for each of the data elements in X,
           all in functional form involving the parameters of the model.  This
           vector is denoted Fvec in the NonlinearLeastSquares class.
 
        4) The display function to be used for plotting the partial and the final
           results if at all such results can be displayed in the form of a 2D or
           3D graphic with Python's matplotlib library.
 
    Finally, the statement in line (D) lets you call on the Levenberg-Marquardt
    algorithm for estimating the model parameters. If you would rather use the more
    basic gradient-descent algorithm for that purpose, your call will look like what
    is shown in line (E).
 
    See the script 
 
                     leven_marq.py
 
    in the ExamplesOptimizedSurfaceFit subdirectory of the distribution for
    further information on what calls you need to sequence together for optimally
    estimating the parameters of an analytically specified surface that you want
    to fit to noisy data.
 
 
USAGE FOR ESTIMATING THE SCENE STRUCTURE WITH A CALIBRATED CAMERA IN MOTION:
 
    You again need to create an instance of the NonlinearLeastSquares class as
    shown in line (F) below.  In addition, you need to create an instance of a
    class like ProjectiveCamera as shown at line (G) below.  If you are writing
    your own code for camera modeling, it would be best if you subclass your code
    from the ProjectiveCamera class in the module.
 
        optimizer =  NonlinearLeastSquares.NonlinearLeastSquares(            #(F)
                                             max_iterations = 400,
                                             delta_for_jacobian = 0.000001,
                     )
        
        camera = ProjectiveCamera.ProjectiveCamera(                          #(G)
                             camera_type = 'projective',
                             alpha_x = 1000.0,
                             alpha_y = 1000.0,
                             x0 = 300.0,
                             y0 = 250.0,
                 )
        camera.initialize()
 
        world_points = camera.make_world_points_random(30)
        world_points_xformed = camera.apply_transformation_to_generic_world_points(world_points, ..... )
 
        ##
        ##  all_pixels = Move the camera to different positions and 
        ##               collect the pixels
        ##
        ##  Now construct the X vector and the corresponding prediction
        ##  vector (Fvec) for nonlinear least-squares:
 
        camera.construct_X_vector( all_pixels )                              #(H)
        camera.construct_parameter_vec_for_calibrated_cameras()              #(I)
        camera.construct_Fvec_for_calibrated_cameras(camera_params_dict)     #(J)
        result = camera.get_scene_structure_from_camera_motion('lm')         #(K)
 
    where the returned result will include the optimal estimations for the scene
    parameters.  Note that the parameter vector constructed in line (I) defines the
    hyperplane over which is defined the cost function.  The goal of nonlinear
    least-squares is return that point on the parameter hyperplane where the cost
    function takes on the least value.  Line (J) constructs the prediction vector in
    terms of the parameters used in line (I). The argument 'lm' in line (K) stands
    for "Levenberg-Marquardt".
 
    See the script 
 
              sfm_with_calibrated_cameras_translations_only.py 
 
    in the ExamplesStructureFromCameraMotion subdirectory of the distribution for
    further information on what calls you need to sequence together for estimating
    the structure of a scene that is being viewed with a moving camera.
 
 
USING SPARSE BUNDLE ADJUSTMENT FOR ESTIMATING THE SCENE STRUCTURE WITH AN UNCALIBRATED CAMERA IN MOTION:
 
    You would again need to construct an instance of the NonlinearLeastSquares and
    an instance of the ProjectiveCamera class as shown in lines (F) and (G) above.
    But now you would need to replace the code in lines (H) through (K) by:
 
        camera.construct_X_vector_for_bundle_adjustment( all_pixels )         #(L)
        camera.construct_parameter_vec_for_uncalibrated_cameras_with_rodrigues_rotations()
                                                                              #(M)
        camera.construct_Fvec_for_bundle_adjustment()                         #(O)
        result = camera.get_scene_structure_from_camera_motion_with_bundle_adjustment()
                                                                              #(P)
 
    The construction of the observation vector X in line (L) is specific to sparse
    bundle adjustment --- in the sense that you first want to list the pixels in all
    the cameras for the first world point, then you want to list the pixels in all
    the cameras for the second world point, and so on.  Ordering the observed data in
    this fashion is necessary to create the sort of block sparsity in the Jacobian
    that is exploited in SBA.
 
    Note that we assume that the intrinsic parameters of the camera are known
    already.  So the goal of the bundle adjustment is to estimate the six extrinsic
    parameters, three for translation of the camera and three (in the form of
    Rodrigues rotations) for the rotation.
        
    See the script 
 
       bundle_adjust_sfm_with_uncalibrated_cameras_translations_only.py
 
    in the ExamplesStructureFromCameraMotion subdirectory of the distribution for
    further information on what calls you need to sequence together for estimating
    both the camera parameters and the scene structure for a scene that is being
    viewed with a moving uncalibrated camera.
 
 
 
INTRODUCTION:
 
    Nonlinear least-squares is a widely used set of algorithms for fitting a model to
    noisy data through the minimization of a cost function.  The observed data is
    represented as a vector (denoted X in this module) and how the model would
    predict each element of X by another vector that is denoted Fvec.  As you would
    expect, each element of Fvec is a function of the parameters of the model.  The
    goal of nonlinear least-squares is to find the best possible values for the
    parameters, best in the sense of minimizing a cost function that is almost always
    the square of the vector norm of the difference between X and Fvec.
 
    The simplest approach to solving a nonlinear least-squares problem is Gradient
    Descent (GD).  The mental imagery that best explains GD is to think of all of the
    model parameters as constituting a hyperplane and the cost function as a surface
    over the hyperplane.  Gradient descent consists of starting at some point in the
    hyperplane, looking straight up at the surface, and walking in the direction of
    the steepest descent on the surface. At each point on the hyperplane, you make
    the size of your next step proportional to the value of the gradient on the
    cost-function surface.  Assuming there are no local minima to trap your progress,
    GD will eventually take you to the point on the hyperplane that is directly below
    the global minimum for the cost function.
 
    Even when getting trapped in a local minimum is not an issue (because, let's say,
    you made a good choice for the starting point), the basic gradient-descent
    algorithm suffers from the shortcoming that the closer you get to the minimum,
    the smaller your steps will be, and, therefore, the slower your progress toward
    the destination.
 
    One way to get round this problem with gradient-descent is to use the
    Gauss-Newton formula to make a direct jump to the minimum --- assuming you are
    already sufficiently close to it.  However, should you inadvertently try to make
    a Gauss-Newton jump when still far from the minimum, your algorithm could become
    numerically unstable and crash.
 
    The algorithm that does a good job of combining the numerical robustness of
    gradient-descent with the speed of Gauss-Newton, while making sure that the
    latter is not invoked too early, is the Levenberg-Marquardt (LM) algorithm. Given
    a start point on the hyperplane spanned by the model parameters, LM starts by
    taking a GD step. In subsequent iterations, before committing itself to a step,
    LM checks whether or not that step can safely be taken with GN.  If not, it
    steers the path toward GD.  However, if the condition for taking the GN step
    safely is satisfied, it goes ahead with that.  In this manner, LM can get to the
    minimum in a small number of steps, often under 10, whereas for the same problem
    and the same start point, GD might take more than a hundred.
    
    The NonlinearLeastSquares class in this module provides implementations for both
    the basic Gradient Descent (GD) algorithm and the more efficient
    Levenberg-Marquardt algorithm for getting to the minimum of a cost function.
    Starting with Version 2.0.0. this class also includes an SBA (Sparse Bundle
    Adjustment) variant of the Levenberg-Marquardt algorithm for optimal computation
    of both the structure and the camera parameters with the data collected by an
    uncalibrated camera in motion.
 
    From a programming standpoint, the most notable feature of NonlinearLeastSquares
    is that it is domain agnostic.  That is, you should be able to use
    NonlinearLeastSquares for solving any problem that requires a GD or an LM
    solution for finding the optimum values for a set of model parameters through the
    minimization of a cost function.
 
    The fact that NonlinearLeastSquares is generic implies that you have to write a
    class of your own for the specific domain in which you wish to use
    NonlinearLeastSquares.  The domain specific class that you create must export to
    NonlinearLeastSquares values for the following options in its constructor:
 
    -- The observed data vector.  This vector is denoted X in the
       NonlinearLeastSquares class.
 
    -- The initial values to be used for the parameters of the model.
 
    -- A vector of the predicted values for each of the data elements in X, all in
       functional form involving the parameters of the model.  This vector is denoted
       Fvec in the NonlinearLeastSquares class.
 
    -- The display function to be used for plotting the partial and the final results
       if at all such results can be displayed in the form of a 2D or 3D graphic with
       Python's matplotlib library.
 
    And, if you wish for NonlinearLeastSquares to use your analytically specified
    partial derivatives in the Jacobian matrix that it needs for the next-step
    calculations, your domain-specific class must also export that matrix to
    NonlinearLeastSquares.  In the absence of a user-supplied Jacobian matrix,
    NonlinearLeastSquares estimates it numerically. [See the implementation code for
    the OptimizedSurfaceFit class supplied with this module --- that class is an
    example of a domain-specific class that uses NonlinearLeastSquares for carrying
    out the needed optimization --- for how your own domain-specific class can
    construct a Jacobian matrix in a functional form and supply it to
    NonlinearLeastSquares.]
 
    The NonlinearLeastSquares class provides several setter methods that your own
    domain-specific class can use to convey the above-mentioned information to
    NonlinearLeastSquares.
 
    To illustrate how you may wish to write your domain specific class, this module
    also comes with two additional classes named OptimizedSurfaceFit and
    ProjectiveCamera, the first for fitting surfaces to noisy data over an xy-plane
    and the second for estimating the structure of a 3D scene from the data recorded
    by a camera in motion.
 
 
METHODS:
 
    Constructing an instance of NonlinearLeastSquares:
 
        optimizer =  NonlinearLeastSquares(   
                         X = None,
                         Fvec = None,
                         num_measurements = None,
                         num_parameters = None,
                         initial_params_dict = None, 
                         initial_param_values_file = None, 
                         jacobian_functionals_array = None,
                         delta_for_jacobian = 0.000001,
                         delta_for_step_size = 0.0001,
                         max_iterations = 200,
                     )
 
        In most usage scenarios, though, you are likely to call NonlinearLeastSquares
        directly with just the following constructor options and let your own domain
        specific class set the other options at run time.  That is, in your own code,
        you will first create an instance of NonlinearLeastSquares through the
        following call to its constructor:
 
        optimizer =  NonlinearLeastSquares(             
                         max_iterations = 400,
                         delta_for_jacobian = 0.000001,
                     )
 
        and subsequently have your own domain-specific class call the various setter
        methods of NonlinearLeastSquares for giving values to its other constructor
        options.  As to which setter methods your own class should call is presented
        in the rest of this section.
 
CONSTRUCTOR OPTIONS:
 
        debug, debug2:
 
            When set to True, it prints out a lot of information at each iteration of
            the nonlinear least-squares algorithm invoked.
 
        display_function:
 
            If the problem you are trying to solve with nonlinear least-squares
            allows for the result of optimization to be visualized, use this
            constructor option to supply a reference to the function you would like
            to be used for such visualization.
    
        Fvec:
        Fvec_BA:
    
            These must be set to a numpy vector (actually a numpy matrix with just
            one column) whose elements are the "predictors" for the corresponding
            observed values in the X and the X_BA vectors, respectively.  (We use the
            symbol 'X' to denote the vector of measured data, as you will soon
            see. X_BA does the same for the bundle-adjustment variant of the basic
            Levenberg-Marquardt algorithm.) Each element of Fvec and Fvec_BA is, in
            general, a function of all of the model parameters.
 
        initial_params_dict:
    
            This is set to a dictionary whose keys are the model parameters and whose
            values the initial values for the model parameters.  The initial values
            for the model parameters specify the point in the parameter hyperplane
            where you want to start the descent to the minimum.
 
        initial_param_values_file:
 
            If your problem involves a very large number of parameters, it may be
            more convenient to place all their initial values in a text file.  Each
            parameter and its initial value must be in a line all by itself.  See the
            file "initial_params_gd2.txt" in the Examples directory for how this file
            needs to be formatted.
 
        jacobian_functionals_array:
 
            If you wish to specify the partial derivatives of the functional elements
            in the Fvec vector with respect to the model parameters, you can supply
            those as a numpy chararray through this constructor option.
 
        num_measurements:
 
            This is simply the number of data values (meaning the number of
            measurements) in X.
    
        num_parameters:
    
            This is the number of model parameters that you wish to calculate with
            nonlinear least-squares.
 
        X:
        X_BA:
 
            This variable must be set to a numpy vector (actually a numpy matrix with
            just one column) whose elements constitute the observed data.  The number
            of elements in X would equal the number of observed data elements that
            you are using for calculating the optimum values of the model parameters.
            Whereas the variable X stores the observations for any regular
            application of the Levenberg-Marquardt algorithm, the variable X_BA
            stores the observed data for the case when you want to use the
            bundle-adjustment variant of Levenberg-Marquardt.  Whereas the generic
            Levenberg-Marquardt places no constraints on how the observations are
            arranged in the vector X, that is not the case with the bundle-adjustment
            variant of the same algorithm.  Both X and M_BA must be a numpy vector,
            meaning a numpy matrix with just one column.
 
 
METHODS:
 
    (1)  grad_descent():
 
         This is the implementation of the basic gradient-descent algorithm.
 
    (2)  leven_marq():
 
         This is the implementation of the Levenberg-Marquardt algorithm mentioned in
         the Introduction section.
 
    (3)  set_debug():
 
         This allow your own domain-specific class to set the 'debug' attribute of an
         instance of NonlinearLeastSquares.
 
    (4)  set_display_function():
 
         Some problem domains allow the result of a nonlinear least-squares
         calculation to be displayed with 2D or 3D graphics.  For example, if the
         goal is to fit an analytical form (in the form of, say, a polynomial) to
         noisy height data over a flat plane and you use the nonlinear least-squares
         algorithm to calculate the best values for the parameters of the analytical
         form, you should be able to visualize the quality of your results by
         displaying both the original noisy data and the model you fit to the data.
         When such a visualization of the results is possible, you can pass the
         definition of the your display function through this setter function.
 
    (5)  set_Fvec():
         set_Fvec_BA()
 
         Each of these methods expects for its main argument a numpy matrix with a
         single column whose each element must be a functional form that predicts the
         corresponding element in the measurement vector X or X_BA.  These functional
         forms will be functions of the model parameters.  The first of the two
         methods is for the regular implementation of the Levenberg-Marquardt
         algorithm and the second for the bundle-adjustment variant of the same.
 
    (6)  set_initial_params():
 
         This method expects for its main argument a dictionary of <key,value> pairs
         in which the keys are the model parameters and the the corresponding values
         the initial values to be given to those parameters.
 
    (7)  set_jacobian_functionals_array():
 
         This method expects for its argument an Nxp matrix of functionals for the
         partial derivatives needed for the Jacobian matrix.  N is the number of data
         elements in the X vector and p is the number of parameters in the model.  To
         elaborate, if you are using nonlinear least-squares to fit an optimal
         surface to noisy height values over the xy-plane, your X vector will be a
         single-column numpy matrix and each row of this vector would correspond to
         one height value at some (x,y) point. The corresponding row in the argument
         jacobian_functionals_array contains p functionals, with each functional
         being a partial derivative of the model functional (with its x and y set
         according to where the height was recorded) with respect to the parameter
         corresponding to the column index.
 
    (8)  set_num_measurements():
 
         This method sets the number of data elements in the X vector in the instance
         of the NonlinearLeastSquares on which the method is invoked.
 
    (9)  set_num_parameters():
 
         This method sets the number of model parameters that will be used in the
         nonlinear least-squares optimization.
 
    (10) set_X():
         set_X_BA():
                
         These two methods set the data vector, in the form of a numpy matrix
         consisting of only one column, in an instance of NonlinearLeastSquares.  The
         first of these for the regular implementation of the Levenberg-Marquardt
         algorithm and the second for the bundle-adjustment version of the same.
         This is the data to which you which you want to fit a given model and you
         want NonlinearLeastSquares to estimate the best values for the parameters of
         the model.
 
OPTIMIZED SURFACE FITTING:
 
    This section presents a class named OptimizedSurfaceFit to illustrate how you can
    use the functionality of NonlinearLeastSquares in your own code.  The goal of
    OptimizedSurfaceFit is to fit a model surface to noisy height data over the
    xy-plane in the xyz-coordinate frame, with the model surface being described by a
    polynomial function.  Here are some examples of such polynomials:
 
           "a*(x-b)**2 + c*(y-d)**2 + e"
 
           "a*x**2 + c*y**2"
 
           "a*x + b*y + c"
    
    where the value returned by the polynomial for given values of the coordinate
    pair (x,y) is the height above the xy-plane at that point.  Given the sort of a
    model surface shown above, the problem becomes one of optimally estimating the
    value for the model parameters from the noisy observed data.  If vector X
    represents the measured data over a set of (x,y) points in the form of a vector
    of observations, we can now write for the cost function:
 
          d^2    =    || X  -  Fvec ||^2
 
    where Fvec is a vector of the predictions as dictated by the model at each of the
    (x,y) point.  That is, each element of the Fvec vector is a prediction for the
    corresponding element of the measurement vector X.  The quantity d^2 is the
    square of the vector norm of the prediction error, meaning the difference between
    the observations in X and the predictions in Fvec.  Given X and Fvec vectors, We
    can call on NonlinearLeastSquares to help us find the best values for the
    parameters of the model surface.
 
    A typical call to OptimizedSurfaceFit's constructor looks like:
 
        surface_fitter = OptimizedSurfaceFit(                         
                                gen_data_synthetically = True,
                                datagen_functional = "7.8*(x - 0.5)**3 + 2.2*(y - 0.5)**2",
                                data_array_size = (16,16),
                                how_much_noise_for_synthetic_data = 0.7,
                                model_functional = "a*(x-b)**3 + c*(y-d)**2",
                                initial_param_values = {'a':2.0, 'b':0.4, 'c':0.8, 'd':0.4},
                                display_needed = True,
                                display_size = (12,8),                 
                                display_position = (500,300),
                                debug = True,
                         )
 
    or, if you wish to also supply the partial derivatives of the model functional
    that can be used by OptimizedSurfaceFit for specifying the the Jacobian matrix to
    NonlinearLeastSquares, like
 
        surface_fitter = OptimizedSurfaceFit(
                                gen_data_synthetically = True,
                                datagen_functional = "7.8*(x - 0.5)**3 + 2.2*(y - 0.5)**2",
                                data_array_size = (16,16), 
                                how_much_noise_for_synthetic_data = 0.5,
                                model_functional = "a*(x-b)**3 + c*(y-d)**2",
                                initial_param_values = {'a':2.0, 'b':0.4, 'c':0.8, 'd':0.4},
                                partials_for_jacobian = {'a':'(x-b)**2', 
                                                         'b':'-2*a*(x-b)', 
                                                         'c':'(y-d)**2', 
                                                         'd':'-2*c*(y-d)'},            
                                display_needed = True,
                                display_size = (12,8),                 
                                display_position = (500,300),
                         )
    
    With regard to the constructor option 'partials_for_jacobian', it is important to
    realize that what is passed to OptimizedSurfaceFit's constructor is not directly
    the Nxp Jacobian matrix (where N is the number of observations in X and p the
    number of parameters in the model).  Instead, it is a set of partial derivatives
    of the model functional with respect to the parameters of the functional.
    However, OptimizedSurfaceFit knows how to translate into an Nxp numpy chararray
    of the functionals needed for the Jacobian.
 
    The constructor options for the OptimizedSurfaceFit class:
 
        data_array_size:                      
 
           The synthetic height that is generated by OptimizedSurfaceFit is over a
           unit square in the xy-plane. Both the x and the y coordinates of this
           square range over the interval 0.0 to 1.0.  If you set this constructor
           option to, say, (16,16), the unit square will be sampled over a 16x16
           grid.
 
        datagen_functional:           
 
            When gen_data_synthetically is set to True, you must supply an algebraic
            expression in the form of a string that OptimizedSurfaceFit can use to
            generate the height data.  Here is an example of what such a string looks
            like:
 
                       "7.8*(x - 0.5)**3 + 2.2*(y - 0.5)**2"
 
            You can call any function inside the string that the Python math library
            knows about.
 
        debug:
 
            This flag is passed on to NonlinearLeastSquares.  When set to True, it
            causes that class to display useful information during each iteration of
            the nonlinear least-squares algorithm.
 
        display_needed:
 
            Fitting optimal surfaces to height data lends itself well to 3D
            visualization.  So if you'd like to see the the surfaces that correspond
            to the optimal values for the model parameters, set this constructor
            option to True.
 
        display_position:
 
            It is set to a tuple of two integers, with the first integer specifying
            the horizontal coordinate and the second the vertical coordinate of the
            upper left corner of the display.  These two coordinate values are with
            respect to the upper left corner of your terminal screen. Horizontal
            coordinates are positive to the right and vertical coordinates positive
            pointing down.  Setting this constructor parameter is optional.  If not
            set, matplotlib will use its default values.
 
        display_size:
 
            It is set to a tuple of two integers, with the first integer specifying
            the width and the second the height of the display.  Setting this
            constructor parameter is optional.  If not set, matplotlib will use its
            default values.
 
        gen_data_synthetically:
 
            If set to True, OptimizedSurfaceFit can generate the measurement height
            data for your synthetically according to the function you specify through
            the constructor option 'datagen_functional'.
 
        how_much_noise_for_synthetic_data:                    
 
           This option controls the amount of noise that is added to the height data
           generated according to the datagen_functional.  The best way to give a
           meaningful value to this construction option is to set to some fraction of
           the largest coefficient in the datagen_functional.
 
        initial_param_values:   
 
           Through this option, you can transmit to OptimizedSurfaceFit your initial
           guesses for the values of the parameters in the model functional.  If you
           cannot think of a guess, you try setting all the parameters to zero.
           OptimizedSurfaceFit conveys your initial values for the parameters to the
           NonlinearLeastSquares class.  You express the initial values in the form
           of a dictionary, whole keys are the name of the parameters and whose
           values the initial values to be given to those parameters.
 
        initial_param_values_file:
 
           If the number of parameters in the problem you are addressing is large, it
           may be more convenient to supply the initial values for the parameters
           through a text file.
 
        measured_data_file:
 
           If the amount of measured data is large, it may be more convenient to feed
           it into the module through a text file.
 
        model_functional:                        
 
           This is the algebraic expression that we want to fit to the noisy height
           data.  OptimizedSurfaceFit will call on NonlinearLeastSquares to estimate
           the best values for the parameters of this algebraic expression.  For
           example, if the model_functional is "a*(x-b)**3 + c*(y-d)**2", the
           NonlinearLeastSquares will find the best possible values for the
           parameters a, b, c, and d --- best in the sense of minimizing the cost
           function described previously.
 
        model_functional_file:
 
           If the model functional is too long and/or too complex to be specified as
           an option directly in a call to the constructor, you can also place it in
           a text file through the model_functional_file option.
 
        optimizer:
 
           Through this constructor option, you can have an instance variable of the
           same name to hold a reference to an instance of NonlinearLeastSquares.
 
        partials_for_jacobian:
 
           Although the NonlinearLeastSquares class can numerically estimate the
           partial derivatives of the element of the Fvec vector with respect to the
           model parameters, with this option you can supply your own analytical
           forms for the partial derivatives that OptimizedSurfaceFit can convert
           into a Jacobian matrix before transmitting it to NonlinearLeastSquares.
 
    Here are the methods defined for OptimizedSurfaceFit:
 
        (1) construct_Fvec():
 
            This method constructs the Fvec vector that NonlinearLeastSquares needs
            for comparing with the measurement vector X. Each element of Fvec is a
            prediction for the corresponding element of X and this prediction is a
            functional form involving model parameters.
 
        (2) construct_jacobian_array_in_functional_form():
 
            This method is used only when the user supplies analytical forms for the
            partial derivatives of the model functional with respect to each of the
            model parameters.  (When the user does not supply such partial
            derivatives, NonlinearLeastSquares estimates the Jacobian through
            numerical approximations.)
 
        (3) display_function()
 
            The problem addressed by OptimizedSurfaceFit lends itself well to
            visualization of the quality of the results returned by
            NonlinearLeastSquares.  With the definition for this method that is
            provided, you can see how well the model parameters estimated by
            NonlinearLeastSquares fit the noisy height data.
 
        (4) gen_data():
 
            This method generates the height data over the xy-plane according to the
            analytical form that is supplied to it as its main argument.  We refer to
            this analytical form as the 'model functional'.
 
        (5) get_initial_params_from_file():
 
            If you want to use model functions that have a large number of
            parameters, it might be easier to place their values in a text file and
            have OptimizedSurfaceFit get it from the file by using this method.
 
        (6) get_measured_data_from_text_file():
 
            If you would like to supply the height data through a text file (rather
            than have the class generate it automatically for you), then this is the
            method to call for reading in the data from the file.  The method assumes
            that that individual data elements are separated by whitespace characters
            (space, tab, newline, etc.).  Since OptimizedSurfaceFit knows about the
            size of the array both along the x-coordinate and along the y-coordinate,
            it knows how to interpret the data in the text file.
 
        (7) get_model_functional_from_file():
 
            If the model functional is too long, you can get OptimizedSurfaceFit to
            read it from a text file by using this option.
 
        (8) set_constructor_options_for_optimizer()
 
            The responsibility of this method is to take all of the user supplied
            information and reconstitute it into a form that is needed by
            NonlinearLeastSquares taking into account the peculiarities of your
            domain.
 
 
ESTIMATING SCENE STRUCTURE WHEN USING A CALIBRATED CAMERA IN MOTION:
 
    This section presents a class named ProjectiveCamera to illustrate how you can
    use the functionality of NonlinearLeastSquares in your own code for estimating
    the structure of a 3D scene from the data recorded by a calibrated camera in
    motion.
 
    To create a simulated structure-from-camera-motion demonstration with this
    module, you must first create an instance the ProjectiveCamera class.  A typical
    call to ProjectiveCamera's constructor looks like:
 
        camera = ProjectiveCamera.ProjectiveCamera(
                             camera_type = 'projective',
                             alpha_x = 1000.0,
                             alpha_y = 1000.0,
                             x0 = 300.0,
                             y0 = 250.0,
                 )
        camera.initialize()
 
    This returns a camera whose optic axis is aligned with the world-Z axis and whose
    image plane is parallel to the world-XY plane. The parameters 'alpha_x' and
    'alpha_y' are for the focal length of the camera in terms of the image sampling
    intervals along the x-axis and along the y-axis, respectively.  The parameters
    'x0' and 'y0' are for the coordinates of the point in the camera image plane
    where the optic axis penetrates the image plane with respect to the origin in the
    image plane (which is usually a corner of the image).
 
        world_points = camera.make_world_points_for_triangle()
        world_points_xformed = camera.apply_transformation_to_generic_world_points( world_points, (0,0,0), (0.0,0.0,5000.0), 1.0)
 
    which generates a triangle defined by its three vertices from a method defined
    for the ProjectiveCamera class and then moves the scene triangle along the optic
    axis of the camera (the world-Z axis) by 5000 units.  After the transformation,
    the three vertices are at the coordinates (3000,3000,5000), (4000,3000,5000), and
    (4000,5000,5000).
 
    Subsequently, you must move the camera to different positions and orientations
    and use the camera matrix constructed by the ProjectiveCamera instance to project
    the world triangle into the camera images. You are going to need the following
    two methods defined for the ProjectiveCamera class for these camera motions:
 
        rotate_previously_initialized_camera_around_x_axis( theta_x_delta )
 
        translate_a_previously_initialized_camera( (0.0,y_motion_delta,0.0) )
 
    At each camera position/orientation achieved with the above two methods, you can
    record the pixels with the following call:
 
        pixels = camera.get_pixels_for_a_sequence_of_world_points( world_points_xformed )
 
    Subsequently, you must make the following call:
 
        construct_X_vector( all_pixels )        
 
    where 'all_pixels' is the set of all the pixel recorded in all the positions of
    the camera.
 
    You would also need to create a Prediction Vector, Fvec, for the observed data
    whose elements are predictor functionals in terms of the scene parameters that
    need to be estimated.  This is achieved with a call like:
 
        construct_Fvec_for_calibrated_cameras( camera_params_dict )
 
    Now you are ready to call the following method:
 
        get_scene_structure_from_camera_motion('lm')
 
    which will invoke the Levenberg-Marquardt method on the NonlinearLeastSquares
    class to estimate the scene structure.
 
    The constructor options for the ProjectiveCamera class:
 
        camera_type:
 
           You can only set this constructor option to 'projective' in the current
           version of the module.  Eventually, I intend to include 'orthographic' as
           another possibility for this option.
 
        alpha_x:
        alpha_y:
 
           These options are for the focal length in terms of the image sampling
           intervals used along the image x-axis and along the image y-axis,
           respectively.
 
        x0:
        y0:
 
           These options are for the coordinates of the point in the camera image
           plane where the optic axis penetrates the image plane with respect to the
           origin in the image plane (which is usually a corner of the image).
 
        camera_rotation:
 
           Using the (roll,pitch,yaw) convention you can specify the rotation for the
           camera in the constructor itself.  However, for experimenting with
           structure-from-camera-motion experiments, it is easier to first construct
           a camera in its generic pose and to then call the rotate and translate
           methods on it in order to move to a different position and orientation.
 
        camera_translation:
 
           Using a triple to indicate displacements along the world-X, world-Y, and
           world-Z, you can specify a translation for the camera in the constructor
           itself.  However, for experimenting with structure-from-camera-motion
           experiments, it is easier to first construct a camera in its generic pose
           and to then call the rotate and translate methods on it in order to move
           to a different position and orientation.
 
 
    Here are the methods defined for ProjectiveCamera:
 
        (1) add_new_camera_to_list_of_cameras():
 
            You will find this utility method useful for enumerating all the
            different camera positions you will be using in a simulated
            structure-from-camera-motion experiment.
 
        (2) apply_transformation_to_generic_world_points()
 
            After you have constructed a scene object (typically just a simple shape
            like a triangle or a tetrahedron), you can call on this method to change
            its position and the pose in the world frame.  The method takes FOUR
            arguments: (1) The scene structure in the form of a list of homogeneous
            coordinates for the world points on the object.  (2) The first is a
            triple that specifies the rotation using the (roll,pitch,yaw)
            convention. (3) The second argument is a triple for the displacement
            along the world-X, world-Y, and world-Z coordinates. (4) The scale factor
            by which you want to expand or shrink the scene object.
 
        (3) construct_Fvec_for_calibrated_cameras(camera_params_dict)
 
            This method constructs the prediction vector Fvec vector that
            NonlinearLeastSquares needs for comparing with the measurement vector
            X. Each element of Fvec is a prediction for the corresponding element of
            X and this prediction is a functional form involving the structure
            parameters.
            
        (4) construct_parameter_vec_for_calibrated_cameras()
 
            This method constructs an ordered list of the SYMBOLIC NAMES to be used
            for each of the coordinates for the scene points that need to be
            estimated.  This list looks like "['X_0', 'Y_0', 'Z_0', 'X_1', 'Y_1',
            'Z_1', 'X_2' ......]"
 
        (5) construct_structure_ground_truth()
 
            This method packages the scene world points in a way that makes it
            convenient to output in your terminal window the estimated coordinates
            for the scene points, the ground-truth value for those coordinates, and
            the initial guesses supplied for them.
 
        (6) construct_X_vector(all_pixels)
 
            As mentioned in the Introduction, we use the notation X to represent a
            vector of all the observed data.  For a structure-from-camera-motion
            problem, the observed data consists of all the pixels in all of the
            camera positions.  This method orders the x- and the y-coordinates of all
            the recorded pixels in the same fashion as the order given to the scene
            points in world-3D.
 
        (7) display_structure()
 
            This method is used to display in the form of a Matplotlib figure the
            following three things simultaneously: the estimated scene structure, the
            actual world points used for the scene object, and the initial guesses
            supplied for those coordinates to the nonlinear least-squares algorithm.
            The three parameters for this method are named
            'structure_points_estimated', 'world_points_xformed', and
            'initial_values_supplied'.
 
        (8) get_all_cameras()
 
            This utility method is convenient for getting hold of all the cameras
            that supplied the data for solving the structure-from- camera-motion
            problem.  We consider an instance of ProjectiveCamera at each of its
            positions in world-3D as a distinct camera.  So if you move the camera
            to, say, 20 different locations, you are in effect using 20 cameras.
 
        (9) get_pixels_for_a_sequence_of_world_points()
 
            For any given camera position, this method applies the corresponding
            camera matrix to each world point, which must be in homogeneous
            coordinates, in the sequence of world points supplied to the method as
            its argument.
 
        (10) get_scene_structure_from_camera_motion('lm')
             get_scene_structure_from_camera_motion_with_bundle_adjustment()
 
            You must call the first of these two methods for estimating the scene
            structure for the calibrated cameras case after you have collected all
            the pixel data from all the different positions of the camera.
            Obviously, before you can call this method, you would need to construct
            the observation vector X from the pixel data the predictor vector Fvec
            from the parameters of the cameras at each of their positions.  And you
            call the second method for doing the same for the case of uncalibrated
            cameras if you want to use the bundle-adjustment version of the
            Levenberg-Marquardt algorithm.
 
        (11) initialize()
    
            This method packs the constructor options supplied to the
            ProjectiveCamera constructor in the form of the camera's intrinsic
            parameter matrix K.  If a translation and/or a rotation is specified for
            the camera through the constructor, those are also incorporated in the
            3x4 camera matrix P put together by this method.
        
        (12) make_world_points_for_triangle()
 
            This method returns a scene object that consists of a triangle in world
            3D. I have found a triangle defined by its three world points to be
            convenient for testing the basic logic of the algorithm for solving a
            structure-from-camera-motion problem.  The triangle returned by this
            method can be subject to any orientation changing and position changing
            transformation.
        
        (13) make_world_points_from_tetrahedron_generic()
 
            Like the previous method, this method returns world points on a
            tetrahedron in world 3D that you can subsequently use for your simulated
            structure-from-moving-camera experiment.
 
        (14) print_camera_matrix()
 
            This utility method is convenient for displaying the 3x4 camera matrix
            for any or all of the positions of the camera.
 
        (15) rotate_previously_initialized_camera_around_world_X_axis()
 
            This method incrementally rotates the camera clockwise around the
            world-X axis by an angle 'theta' in degrees that is supplied to the
            method as its argument.
 
        (16) rotate_previously_initialized_camera_around_world_Y_axis()
 
            This method incrementally rotates the camera clockwise around the world-Y
            axis by an angle 'theta' in degrees that is supplied to the method as its
            argument.
 
        (17) set_constructor_options_for_optimizer( optimizer )
             set_constructor_options_for_optimizer_BA( optimizer )
 
            A ProjectiveCamera instance uses these method to pass on to
            NonlinearLeastSquares all the information needed by the latter (such as
            the observed data vector X or X_BA and the prediction vector Fvec or
            Fvec_BA) for constructing an optimum estimate of the scene structure.
            The argument 'optimizer' that this method takes is an instance of
            NonlinearLeastSquares.
 
        (18) set_initial_values_for_params()
 
            Every nonlinear least-squares algorithm needs a starting guess for
            whatever it is that is being estimated.  In most cases, you would
            construct a random guess for the parameters and supply those values to
            this method in the form of a dictionary in which each key is the symbolic
            name of one of the parameters being estimated and the value a random
            guess for the parameter.
 
        (19) set_params_list( params_arranged_list )
            
            You will use this method to pass on to the instance of ProjectiveCamera
            an ordered list of the parameters you want estimated with nonlinear
            least-squares.
 
        (20) translate_a_previously_initialized_camera()
 
            This method incrementally displaces the camera by 'translation' that is
            supplied to it as its argument. The argument 'translation' consists of a
            triple of real numbers that stand for a displacement along the world-X,
            along the world-Y, and along the world-Z.
 
        (21) displayImage(self, argimage, file_name_for_save="")
 
            This function, added in Version 2.0.2, is for constructing a camera image
            of an object defined virtually in the world XYZ frame.  I have used this
            function in the script flyby_photos_of_3D_planar_shape.py in the
            ExamplesStructureFromCameraMotion directory for constructing flyby
            sequence of images of a planar scene in the YZ-plane of the world frame.
 
 
ESTIMATING SCENE STRUCTURE AND CAMERA PARAMETERS WHEN USING AN UNCALIBRATED CAMERA IN MOTION:
    
    The problem of scene construction becomes a lot more challenging when the camera
    in motion is uncalibrated.  When I say uncalibrated, I mean uncalibrated with
    respect to its extrinsic parameters.  We assume in all cases in this document
    that the intrinsic parameters of the camera are known.
 
    What makes structure estimation more complicated in this case is that each new
    camera position adds 6 additional variables to the overall parameter space, 3 for
    translation and 3 for rotation as expressed by the Rodrigues parameters.  Let's
    say you have M camera positions and N structure points.  That would call for (6*M
    + 3*N) parameters to be estimated by the Levenberg-Marquardt algorithm because
    each camera has six external pose parameters associated with it and each
    structure point is defined by its three world coordinates.  Assuming that you can
    see all of the structure points in all the cameras, all of the pixels recorded in
    all of the camera position would constitute a 2*N*M dimensional observation
    vector X.  Therefore, your Jacobian will be of size [(2*N*M) x (6*M+3*N)].  For
    an example, say you have 20 camera positions and 100 structure points, your
    Jacobian will be of size 4000x420.  Calculating a Jacobian of this size and
    multiplying it by its transpose could take many more than a few minutes on a run
    of the mill machine.  And having to calculate the Jacobian multiple times in the
    iterative framework of nonlinear least-squares estimation could test your
    patience as you are waiting for the results.
 
    Fortunately, this otherwise long execution time can be significantly shortened if
    you take advantage of the sparsity of the Jacobian when using uncalibrated
    cameras.  As to why the Jacobian is sparse, consider the fact that each row of
    the Jacobian is a partial derivative of the prediction for one observation with
    respect to ALL the parameters.  Obviously, when you take the partial derivative
    of the prediction function for a pixel in one specific camera with respect to the
    parameters of all other cameras, all those partial derivatives will be zero.  The
    implementations of the Levenberg-Marquardt algorithm that take advantage of the
    sparsity of the Jacobian are commonly referred to as "Sparse Bundle Adjustment"
    or just "Bundle Adjustment".
 
    This logic that this module uses for exploiting the sparsity of the Jacobian is
    based on the paper "SBA: A Software Package for Generic Sparse Bundle Adjustment"
    by Manolis Lourakis and Antonis Argyros that appeared in ACM Transactions on
    Mathematical Software, March 2009.
    
    In honor of these two authors, I have named all of the block Jacobian submatrices
    that take partial derivatives of the predictors with respect to the camera
    parameters as the Argyros matrices.  And I have named all of the block Jacobian
    matrices that do the same with respect to the structure variables as Lourakis
    matrices.
 
    In this module, the "bundle-adjustment" version of Levenberg-Marquardt algorithm
    can be found in the method
 
           bundle_adjust()
 
    of the NonlinearLeastSquares class.  As you will notice, the observation vector
    that I called X in the implementation of the method leven_marq() is now called
    X_BA in the implementation of the method bundle_adjust().  The reason has to do
    with the fact that the basic Levenberg-Marquardt algorithm that you see in
    leven_marq() does not care how the observed data is arranged in the vector
    X. Typically, in my personal computer vision code for multi-camera situations, I
    arrange the data by the camera.  That is, I group together all of the pixels
    recorded in one camera, followed by all of the pixels in the second camera, and
    so on.  That is how the observed data (and, therefore, also the prediction vector
    Fvec) is arranged in the method leven_marq().  Unfortunately, that ordering of
    the data does not work for sparse bundle adjustment.  To maximally exploit sparse
    bundle adjustment, you must order the observation vector so that your first group
    together the pixels for what we may refer to as the first structure point, to be
    followed by the pixels in all the cameras for the second structure point, and so
    on.  Because this ordering is critical to the SBA algorithm described in the
    paper by Lourakis and Argyros, I have denoted the observation vector X_BA in the
    method bundle_adjust().  The order in which the predictor functionals are placed
    in the predictor vector Fvec must correspond to the order used in X_BA.  So the
    prediction vector in bundle_adjust() is named Fvec_BA.
 
    About how to invoke the bundle-adjustment variant of Levenberg-Marquardt in your
    own code, note the following methods of the ProjectiveCamera class that are
    specially meant for that purpose:
 
        (1)  construct_X_vector_for_bundle_adjustment()
 
             It is this method's job to arrange the observed data (meaning the pixels
             in the camera images) in the special order that is needed by the
             bundle-adjustment algorithm.
 
        (2)  construct_parameter_vec_for_uncalibrated_cameras_with_rodrigues_rotations()    
 
             When using uncalibrated cameras, the parameter vector must include three
             translational and three rotational parameters for each camera.  These
             are in addition to the three (X,Y,Z) parameters for each structure point
             being tracked in the cameras.  The job of this method is to synthesize
             this parameter vector.
             
        (3)  construct_Fvec_for_bundle_adjustment()
 
             The ordering used for the observed data (meaning the pixels in the
             cameras) in the X_BA vector must also be used in the prediction vector
             when you use the bundle-adjustment variant of the basic
             Levenberg-Marquardt algorithm.  This method creates the needed
             prediction vector.
       
        (4)  get_scene_structure_from_camera_motion_with_bundle_adjustment()
 
             It is this method of the ProjectiveCamera class that invokes the
             bundle_adjust() method of the NonlinearLeastSquares class.
 
 
THE ExamplesOptimizedSurfaceFit DIRECTORY:
 
    See the 'ExamplesOptimizedSurfaceFit' directory in the distribution for examples
    of how you can use the NonlinearLeastSquares class for solving optimization
    problems.  These examples are based on the domain specific class
    OptimizedSurfaceFit that knows about fitting model surfaces to noisy height data
    over a flat plane.  You will see the following four scripts in this directory:
 
        leven_marq.py
 
        grad_descent.py    
 
        leven_marq_with_partial_derivatives.py
 
        grad_descent_with_partial_derivatives.py
 
    For the first two scripts, the NonlinearLeastSquares instance used will estimate
    the needed Jacobian matrix through appropriate numerical approximation formulas
    applied to the elements of the Fvec vector.  On the other hand, for the third and
    the fourth scripts, your own domain-specific class must construct the Jacobian
    matrix, in the form of an array of functions. In the case of the domain-specific
    class OptimizedSurfaceFit that comes with this module, this Jacobian matrix is
    constructed from the user-supplied partial derivatives for the model functional.
 
    In order to become familiar with the NonlinearLeastSquares class, you might wish
    to play with the four scripts listed above by:
 
    -- Trying different functional forms for the 'datagen_functional' for
       different shaped surfaces.
 
       When you change the algebraic form of 'datagen_functional' for the
       OptimizedSurfaceFit class, make sure that you also change the algebraic form
       supplied for 'model_functional'.  Note that nonlinear least-squares can only
       calculate the parameters of a model functional that best fit the noisy height
       data; it cannot conjure up a new mathematical form for the surface.  So the
       basic mathematical form of the 'model_functional' must be the same as that of
       the 'datagen_functional'.
 
    -- Trying different degrees of noise.  
 
       As mentioned elsewhere, when you supply a numerical value for the constructor
       option 'how_much_noise_for_synthetic_data' for the OptimizedSurfaceFit class,
       the number you enter should be in proportion to the largest numerical
       coefficient in the 'datagen' functional.  Change this numerical value and see
       what happens to the quality of the final results.
 
    -- Try different values for the initial values of the model parameters.
 
       Since, depending on where the search for the optimum solution is started, all
       nonlinear least-squares methods can get trapped in a local minimum, see what
       happens when you change these initial values.
      
    -- Try different algebraic expressions for the 'model_functional' constructor
       option for the OptimizedSurfaceFit class.  But note that if you change the
       algebraic form of this functional, you must also change the algebraic form of
       the 'datagen_functional' option.
 
    -- Try running the example with and without the partial derivatives that are
       supplied through the 'partials_for_jacobian' option for the
       OptimizedSurfaceFit class.
 
 
THE ExamplesStructureFromCameraMotion DIRECTORY:
 
    See the 'ExamplesStructureFromCameraMotion' directory in the distribution for the
    following three example scripts that show how you can use the
    NonlinearLeastSquares module for estimating the structure of a 3D scene from the
    images recorded by a moving camera:
 
        sfm_with_calibrated_cameras_translations_only.py
 
        sfm_with_uncalibrated_cameras_translations_only.py
 
        bundle_adjust_sfm_with_uncalibrated_cameras_translations_only.py
 
    The first script listed above is for estimating the scene structure with a
    calibrated camera in motion.  As you play with this method, make sure you change
    the level of noise in the initial values supplied for the structure parameters to
    be estimated.  As you will see, the method works even when the initial values for
    the parameters are far from their true values.  Note that the ProjectiveCamera
    class makes it easy to specify calibrated cameras.  The constructor of the class
    first gives you a camera for which you can specify the internal and the external
    parameters through the constructor options. Subsequently, you can apply
    translational and rotational transformations to the camera to move it to
    different locations in world 3D.  Since the 3x4 camera matrices for all these
    positions of the camera are known, you end up with a set of fully calibrated
    cameras for experimenting with structure-from-motion simulations.
 
    The second and the third scripts listed above are for the case of uncalibrated
    cameras, with the former a straightforward application of the Levenberg-Marquardt
    algorithm and the latter a bundle-adjustment variant of the same.  Logically,
    both these methods must return identical answers.  (If you encounter a case when
    the two do not return the same answer, please send a bug report to me. I'd
    appreciate that very much.)
 
    Just to give you an idea of the speed-up you will get with bundle-adjustment,
    when I run the second script listed above on my laptop, it takes about 15 minutes
    for the number of structure points and the number of camera positions used in
    that script.  For exactly the same number of structure points and the camera
    positions, the third script takes only a couple of minutes.  You can only imagine
    the speed-up you will get with a C-based library for sparse bundle adjustment ---
    such as the "sba" library mentioned in the paper by Lourakis and Argyros that I
    cited earlier in this documentation page.
 
    Starting with Version 2.0.2, an additional script in this directory is
 
        flyby_photos_of_3D_planar_shape.py
 
    that allows you to create a flyby sequence of images by "flying" a virtual camera
    (created as an instance of the ProjectiveCamera class) over a scene in the world
    XYZ frame.  The specific example in the above script is for a planar shape in the
    world YZ-plane.
 
 
CAVEAT
 
    Note the bundle-adjustment variant of the Levenberg-Marquardt algorithm that you
    see in the bundle-adjust() method of the NonlinearLeastSquares module is meant
    for just educational purposes.  Being pure Python, it cannot compete with highly
    optimized C-based code you will see in, say, the "sba" library mentioned in the
    article by Lourakis and Argyros that I have cited earlier in this documentation.
    So if your needs for nonlinear least-squares for estimating the structure and the
    camera parameters are primarily of the production variety, go directly to those
    publicly available libraries.
 
 
INSTALLATION:
 
    The NonlinearLeastSquares class was packaged using setuptools.  For installation,
    execute the following command-line in the source directory (this is the directory
    that contains the setup.py file after you have downloaded and uncompressed the
    package):
 
            sudo python setup.py install
    and/or
            sudo python3 setup.py install
 
    On Linux distributions, this will install the module file at a location that
    looks like
 
             /usr/local/lib/python3.6/dist-packages/
 
    If you do not have root access, you have the option of working directly off the
    directory in which you downloaded the software by simply placing the following
    statements at the top of your scripts that use the NonlinearLeastSquares class:
 
            import sys
            sys.path.append( "pathname_to_NonlinearLeastSquares_directory" )
 
    To uninstall the module, simply delete the source directory, locate where the
    NonlinearLeastSquares module was installed with "locate NonlinearLeastSquares"
    and delete those files.  As mentioned above, the full pathname to the installed
    version is likely to look like
    /usr/local/lib/python2.7/dist-packages/NonlinearLeastSquares*
 
    If you want to carry out a non-standard install of the NonlinearLeastSquares
    module, look up the on-line information on Disutils by pointing your browser to
 
              http://docs.python.org/dist/dist.html
 
 
BUGS:
 
    Please notify the author if you encounter any bugs.  When sending email, please
    place the string 'NonlinearLeastSquares' in the subject line.
 
 
ACKNOWLEDGMENTS:
 
    Version 2.0.3 was a result of Somosmita Mitra discovering a couple of bugs in my
    Levenberg-Marquardt code for Nonlinear Least Squares. She not only discovered the
    bugs but also supplied the fixes for them.  The bugs were in the updating of the
    residual error and the damping coefficient from one iteration to the next of the
    algorithm.  Somosmita is finishing her PhD in Purdue RVL and, at the moment, is
    the TA for the Computer Vision class that I teach.  Thanks, Somosmita!
 
 
ABOUT THE AUTHOR:
 
    Not too long ago, the author, Avinash Kak, finished a 17-year long "Objects
    Trilogy Project" with the publication of the last book "Designing with Objects"
    by John-Wiley. If interested, visit his web page at Purdue to find out what this
    project was all about. You might like "Designing with Objects" especially if you
    enjoyed reading Harry Potter as a kid (or even as an adult, for that matter).
 
    For all issues related to this module, contact the author at kak@purdue.edu
 
    If you send email, please place the string "NonlinearLeastSquares" in your
    subject line to get past the author's spam filter.
 
 
 
COPYRIGHT:
 
    Python Software Foundation License
 
    Copyright 2026 Avinash Kak
 
@endofdocs

Imported Modules
		glob itertools math numpy numpy os re scipy sys

Classes
		builtins.object NonlinearLeastSquares   class NonlinearLeastSquares(builtins.object)     NonlinearLeastSquares(*args, **kwargs)       Methods defined here: __init__(self, *args, **kwargs) constructor bundle_adjust(self, *args, **kwargs) This is an implementation of the "bundle adjustment" version of the Levenberg-Marquardt algorithm for nonlinear least-squares.  Bundle adjustment takes advantage of the sparsity of the Jacobian that one sees in  applications such as estimating the scene structure with the images recorded with uncalibrated cameras.   The implementation shown here is based on the now celebrated paper "SBA: A Software Package for Generic Sparse Bundle Adjustment" by Manolis Lourakis and Antonis Argyros that appeared in ACM Transactions on  Mathematical Software, March 2009.   This function is used in the following scripts        bundle_adjust_sfm_with_uncalibrated_cameras_translations_only.py    in the ExamplesStructureFromCameraMotion directory of the distribution.     Note that since bundle_adjust() is written in a generic manner, it is not called directly by the example scripts listed above.  As shown in the scripts, you need to first construct an instance of the class ProjectiveCamera that is a co-class of the main module class NonlinearLeastSquares in the distribution. grad_descent(self) This is an implementation of the basic gradient descent algorithm for nonlinear least-squares as described  in my Lecture 13 notes at the lecture-notes website for Purdue University's ECE661: Computer Vision leven_marq(self) This is an implementation of the Levenberg-Marquardt algorithm for nonlinear least-squares as described  in my Lecture 13 notes at the lecture-notes website for Purdue University's ECE661: Computer Vision   This function is called by the following scripts in the ExamplesOptimizedSurfaceFit directory of the distro:               leven_marq.py               leven_marq_with_partial_derivatives.py   It is important to note that the above two scripts do NOT call leven_marq() function directly.  AS  stated in the main document page for the NonlinearLeastSquares module, the code in the file  NonlinearLeastSquares.py is written in a domain agnostic manner.  So you need a domain adaptation  class that knows how to package the arguments for calling leven_marq().  For the case of the two scripts listed above, that domain-specific class in the distro is OptimizedSurfaceFit.  It is a co-class of  the main module class NonlinearLeastSquares in the distribution.   The function leven_marq() defined here is ALSO used in the following two scripts       sfm_with_calibrated_cameras_translations_only.py     sfm_with_uncalibrated_cameras_translations_only.py   in the ExamplesStructureFromCameraMotion directory of the distribution.  Again, these example scripts do NOT directly call the leven_marq() function.  As shown in the two scripts, they must first construct an instance of the ProjectiveCamera class that knows how to bundle the arguments together for calling the leven_marq() function. set_Fvec(self, Fvector) This method supplies the NonlinearLeastSquares class with the prediction vector whose each element is a functional form of the prediction in the observed data vector X.  Note that  Fvec is a column vector --- meaning a numpy matrix with just one column. set_Fvec_BA(self, Fvector_BA) You need to call this method for providing the NonlinearLeastSquares class with the prediction vector if you are going to be using the bundle-adjustment capabilities of the class. set_X(self, X) set_X_BA(self, X_BA) set_debug(self, debug) set_display_function(self, display_function) set_initial_params(self, initial_params_dict) set_jacobian_functionals_array(self, jacobian_functionals_array) This method expects for its argument an Nxp matrix of functionals for the partial  derivatives needed for the Jacobian matrix.  N is the number of measurements in the X vector and p is the number of parameters in the model.  If you are using nonlinear least-squares to fit optimal surfaces to noisy measurements over the xy-plane, each element of the X vector would correspond to one such measurement at some (x,y) coordinates. And an element the argument jacobian_functionals_array chararray would correspond to the partial derivative of the model functional that already has incorporated the (x,y) coordinates corresponding to that row and that is  a partial derivative of the model with respect to the parameter corresponding to the column. set_num_measurements(self, how_many_measurements) set_num_parameters(self, how_many_parameters) set_params_arranged_list(self, params_list) set_params_ordered_list(self, params_list) set_problem(self, prob) If you are using this module to solve structure from camera motion (sfm) problems, use this method to set 'self.problem' to 'sfm_N' where 'N' is the number of world points you are tracking.   This  is needed because sfm needs the specialized display function defined for the ProjectiveCamera class. Data descriptors defined here: __dict__ dictionary for instance variables (if defined) __weakref__ list of weak references to the object (if defined)

		__author__ = 'Avinash Kak (kak@purdue.edu)' __copyright__ = '(C) 2026 Avinash Kak. Python Software Foundation.' __date__ = '2026-October-12' __url__ = 'https://engineering.purdue.edu/kak/distNonlinearLeastSquares/NonlinearLeastSquares-2.0.4.html' __version__ = '2.0.4'

Author
		Avinash Kak (kak@purdue.edu)