Programming tips for implementation of a non-linear optimization solver (Part-I: Linear algebra)
I. Introduction
In this article and in the series of upcoming articles, I will explain you about a few programming tips that I have been integrating for a custom non-linear optimization solver that I am developing as a pet project. Before getting into the thick of things, let me first explain a few basics on optimization and why I chose to develop my own solver (which I believe is a fair reason to explain considering the fact that a lot of pre-existing solvers exists in the market).
Optimization is a branch of applied mathematics which deals with computing the optimal solutions for a problem which are subjected to a few requirements and/or resource constraints. In short, the idea is to either analytically or algorithmically obtain the best solution for some required objectives over a pool of all other possible solutions. This idea of computing the optimal solution(s) for an optimization problem underpins most of the design, decision, prediction and estimation problems that pops up in almost all domains of engineering problems. Examples include a) computing the best possible control command in control engineering, b) tuning the optimal parameters for a neural network, c) finding the best parameters for mechanical/electrical design problems etc. and so on and so forth. Every single engineering device that we see around us has some element of optimization in it, either this is enforced in the design phase or execution phase. This ubiquitous and applicability nature of the field caught my attention and thereby I embarked on a journey to implement my own optimization solver. Of course, I have just begun this journey, however in due course of my coding I realized a few programming tips which I would like to share with you that might be useful for your application.
Despite there are several nonlinear optimization solvers (Open source or proprietary) such as CPLEX, Gurobi, IPOPT, KNITRO etc. a common trobe that one might lob at me would be “Why re-invent the wheel?”. It is true that I am re-inventing the wheel, but at the same time I am also mastering the skill of “how to build a wheel from scratch”, which I consider far more precious than utilizing pre-existing tools. With this in mind, I have made it as a strict policy to never read the source codes of existing solvers which might hamper my creative potential and style of coding.
The next and final issue which I would like to talk about is the platform which I ought to work on. The simplest choice which I could have settled for would be Python. However, one of my goals was also to focus on high performance and fast execution of the solver, which is the most coveted feature for real-time embedded computations. With this in mind, the natural choice would be either C or C++ programming language. C language is undoubtedly the best language when it comes to execution speed, however without features such as classes, templates, metaprogramming etc., this task might get time-consuming due to the fact that a few codes have to be reimplemented e.g. for different datatypes, a new function have to be written. With templates, one could obviate this trouble. So weighing all pros and cons between C and C++, I decided to go with C++. With Modern C++14/17/20, there are several features that one could utilize to code efficiently and effectively.
When coding a nonlinear optimization solver, I realized that there are so many tools that I have to build on my own such as a Linear Algebra engine, Symbolic engine, Numerics etc. Therefore, I might address the several coding tips in a nonlinear fashion and hopefully not place things in a disarray.
II. Problem formulation
There are several types of optimization problem formulations based on the nature of the decision variables or optimization variables, objectives and constraints functions. Some of the common types of distinctions in optimization problem formulations are
- Continuous vs discrete optimization problems
- Convex vs Non-convex optimization problems
- Finite vs Infinite optimization problems
- Derivative based vs derivative-free optimization problems
For more details on these optimization problems refer [2][3]. Non-linear optimization problem could be simply termed as a union of Non-convex, derivative based and continuous optimization problems. A simple nonlinear optimization problem is given as
where f,g,h are objective, equality constraints and inequality constraints of the nonlinear optimization problem and x denotes the decision variables (shown in Fig2.). It is important to note that the objective and the constraints functions are continuous and twice differentiable and also the decision variables are in the realm of continuous spaces. There are several resources for learning the theoretical properties such as stability, convergence and practical implementation of a nonlinear optimization problem, however, in this article I have chosen the following the books [2][3].
III. Linear algebra
With the basic nonlinear optimization problem formulation as shown in the previous section, there are several numerical methods available for solving it such as Interior point methods (IPM), Sequential quadratic programming (SQP) etc. However, irrespective of the aforementioned methods, for applying a derivative based method, it is important to compute the Jacobians (first order derivatives) or Hessians (second order derivatives) which could be either a vector or a matrix. It is important to perform certain basic operations on it such as addition, multiplication, dot product etc with these matrices. Therefore in this section, as an example, I would like to discuss two different C++ implementations of Matrix-Matrix addition i.e. a naïve approach and a smart approach (using Expression templates). This example could be further extended to Matrix-Vector, Vector-Vector operations etc. which I leave it to the reader to implement it.
a) Coding a Matrix-Matrix addition (Naïve approach)
The code for the naive matrix addition is shown in Fig3. (NaiveMatrix.h). The implementation is a templatized MAT class, where the matrix-matrix addition is implemented with overloading the “+” operator.
A closer look into the operator “+” overloading function is shown in Fig4. The implementation is straight forward. A result matrix of appropriate dimensions is created and in two forloops the corresponding elements are added and stored in the result matrix at appropriate indices. Finally the result matrix is returned. Simple!
a) Coding a Matrix-Matrix addition (Smart approach)
The code for the smart matrix addition is shown in Fig4. (SmartMatrix.h (1st part)). Now, for people who are not comfortable or familiar with template metaprogramming, I know this might not be very intuitive. However, I will try to explain it. Firstly, I implemented an interface namely IMatrix which implements a few methods as show below. The IMatrix is actually a Curiously Recurring Template Pattern (CRTP) interface and the MAT class is a templatized version which inherits IMatrix with the type Matrix<T>. The body of the Matrix is virtually the same as the MAT class shown in Fig3.
The MatrixSum class is shown in Fig5. where the operator “+” is overloaded which returns a MatrixSum class. The MatrixSum class follows the CRTP paradigm just as the Matrix class shown in Fig4.
I know this is not a simple implementation, its convoluted!. So let me try to unpretzel this in simple terms.
How the naive approach works?
In the naive example say if you want to add three matrices D = A + B + C, what really happens is the compiler adds (A+B) and stores in a temporary variable (MAT<T>) and then this temporary variable gets added with C. Once again a temporary variable is created and the values are stored into the D matrix.
How the smart approach works?
In case of the smart approach, when the three matrices are added i.e D = A + B + C, an object of the following template structure is returned without creation of any temporary variable.
MatrixSum<MatrixSum<Matrix<double>, Matrix<double>>, Matrix<double>>
When forcing the evaluation of the generated object, the operator()(size_t, size_t) is invoked, which further invokes the nested operator()(size_t, size_t) inside the object’s template structure and this is repeated recursively and finally the addition is completed. This method is also known as expression templates. For more details refer [4][5].
IV. Results
The above two approaches were implemented on Raspberry Pi 4 Model B board for addition of three matrices. The matrix dimensions were 1000x1000, i.e. a million elements and respective execution time for both the approaches were recorded. The matrices were randomly populated each time and executed for 9 iterations. The following table shows the execution times
The average execution time for the naive and smart approaches are 46715.222 microseconds and 21857.78 microseconds. There is a 2.13x better performance for the smart approach compared to the naive approach. However, with a larger expression than addition of three matrices a well pronounced difference could be obtained.
V. Conclusion
In this article, I have demostrated that the naive approach for matrix addition pales in comparison with respect to the smart approach which employs the expression templates pattern. For non-linear optimization problems, linear algebra is an essential tool and matrix operations as mentioned in this article forms the building blocks. Therefore by optimizing the code, optimal solutions could be computed at a much faster pace.
References
[1] Amini, Alexander, et al. “Spatial uncertainty sampling for end-to-end control.” arXiv preprint arXiv:1805.04829 (2018).
[2] Nocedal, Jorge, and Stephen J. Wright, eds. Numerical optimization. New York, NY: Springer New York, 1999.
[3] Bazaraa, Mokhtar S., Hanif D. Sherali, and Chitharanjan M. Shetty. Nonlinear programming: theory and algorithms. John Wiley & Sons, 2013.
[4] Vandevoorde, David, and Nicolai M. Josuttis. C++ Templates: The Complete Guide, Portable Documents. Addison-Wesley Professional, 2002.