As part of a placement earlier this year, I was given two sets of functions to play with. The first set was implemented in both MATLAB and FORTRAN and formed a root-finding algorithm. The second was implemented in only MATLAB but promised shorter running times and greater accuracy. The task was to compare the two and determine whether implementing the new methods in FORTRAN would be worth the effort.
I used an existing suite made up of scenarios involving large polynomials to test the functions. Unfortunately, the order in which the expected roots were stored was not necessarily the same as those that were produced from the root finders. As the roots were not guaranteed to be real numbers, this could not be solved by simply placing them in order of magnitude. To show why this is the case, we can use the Argand plane which represents complex numbers as coordinate points. The real part is plotted along the x-axis and the complex part is plotted against the y-axis, e.g. 1 + i is plotted at (1, 1) and 3-2i is plotted at (3, -2).
Now we see that for any given magnitude, an infinite number of complex numbers lie on the circumference of a circle with the same radius in the Argand plane. For example, numbers such as 1+i and -1-i would be considered equivalent – clearly undesired behaviour!
A more involved method of matching the computed roots to their closest expected root was needed.
We can use the Argand plane to represent the expected and outputted roots as sets of points. In this representation, we have a matching problem in which we wish to minimise the distance between each of the outputted roots and their nearest expected neighbour.
To do this, I chose the Kuhn-Munkres algorithm. Although it’s of O(n^3) complexity, it was fast enough for the small numbers of roots considered (orders of 100) and was quick to implement.
The first step was to build up a matrix of values to iterate over. This was constructed by setting columns to represent expected values and rows to represent the outputted values. I calculated the distance (in the Argand plane) between each outputted root and every expected root and inserted the distances into the corresponding positions in the matrix. The Kuhn-Munkres algorithm was then applied to this matrix to give the perfect matching of outputted and expected roots and the distances between them.
These distances could then be used to compare the accuracy of the solutions. Although there are many ways this could have been done, I chose simply to use the mean distance and the largest distance. These were compared between each implementation in order to determine which was better. The results showed that the new implementation was both more accurate and faster across the majority of tests and so the re-write went ahead!
- At first glance, tasks may seem more simple than they really are. Nothing should be assumed before gaining a full understanding of the problem domain.
- Unit tests change with the code base – and that’s okay. An alternative solution would have been to order the roots during calculation. Ordering roots in the unit tests allowed the production code to stay slim and without unnecessary operations.