Iterative algorithms – Convergence rates

In optimization and any other iterative numerical algorithms, we are interested in having convergence estimates for all algorithms. We are not only interested in showing that the error goes to , although that is essential, but also at what speed the error goes to . An algorithm which converges quicker is a better choice in all cases!

Generally, there are two points of view for convergence: convergence in terms of where is the number of iterations or having a comparison between the error at iteration and the error at the iteration . I will take the latter point of view below.

To fix the ideas, denote the error estimate for , where is the current iterate and is the solution to the problem.

We have the following standard classification:

• linear convergence: there exists such that
  the constant is called the convergence ratio  it is easy to show that , so in particular .
• sublinear convergence:  but is not linearly converging
• superlinear convergence:  with any positive convergence ratio
   sufficient condition:
• {convergence of order }: there exists such that for large enough
  is called the order of convergence 
  the case has a special name: quadratic convergence

To underline the significant difference between linear and superlinear convergence consider the following examples: Let . Then:

• converges linearly to zero, but not superlinearly
• converges superlinearly to zero, but not quadratically
• converges to zero quadratically

Quadratic convergence is much faster than linear convergence.

Among optimization algorithm, the simpler ones are usually linearly convergent (bracketing algorithms: trisection, Golden search, Bisection algorithm, gradient descent). Algorithms involving higher order information or approximation are generally superlinear (Secant method, Newton, BFGS or LBFGS in higher dimensions etc.).

There is a huge difference between linear convergence and super-linear convergence. If a faster algorithm is available using it is surely useful!

Golden search algorithm – efficient 1D gradient free optimization

Bracketing algorithms for minimizing one dimensional unimodal functions have the form:

• Suppose is an interval containing the minimizer ;
• Pick in and evaluate the function at these points.
• If then choose
• If then choose
• Stop the process when has a length smaller than a given tolerance or a given maximum number of function evaluations is reached.

The simplest algorithm corresponds to choosing which divide the interval into three equal parts. This algorithm can be called trisection algorithm. Below you can see the intervals and intermediary points for a few iterations in the trisection algorithm.

Bracketing intervals and intermediary points for trisection algorithm

Optimizing a 1D function – trisection algorithm

Optimization problems take the classical form

Not all such problems have explicit solution, therefore numerical algorithms may help approximate potential solutions.

Numerical algorithms generally produce a sequence which approximates the minimizer. Information regarding function values and its derivatives are used to generate such an approximation.

The easiest context is one dimensional optimization. The basic intuition regarding optimization algorithms starts by understanding the 1D case. Not all problems are easy to handle for a numerical optimization algorithm. Take a look at the picture below:

photo from Ziv Bar-Joseph

Find coefficients of trigonometric polynomial given samples

Suppose the values of some trigonometric polynomial

are known for a set of distinct angles , . 

Question: Recover the coefficients of the trigonometric polynomial that verify when or which best fits the values in the sense of least squares when . 

Answer: Obviously, since there are unknowns, we need at least equations, therefore we restrict ourselves to the case . The equalities produce a set of equations which has at most a solution when , provided the rank of the matrix of the system is . Define the function

This function is a sum of squares, which has zero as a lower bound. The function can be written using norms as where

A straightforward computation shows that the gradient of is

The matrix is invertible, provided all angles are distinct. This is a direct application of the formula of the Vandermonde determinant: using operations on columns you can recover the Vandermonde matrix corresponding to . Therefore, one can always solve the system when and will minimize . In particular, when the minimum will be equal to zero and the coefficients of the trigonometric polynomial verifying will be found. When the best fit, in the sense of least squares, of the values with a trigonometric polynomial will be found.

Below, you can find a Python code which solves the problem in some particular case.

import numpy as np
import matplotlib.pyplot as plt

N = 5 # coeffs
M = 2*N+1 # M>=2N+1 samples

# function to be approximated
def fun(x):
    return np.sin(x+np.sqrt(2))+0.3*np.sin(5*x-np.sqrt(3))+0.1*np.sin(8*x-np.sqrt(7))

thetas =np.linspace(0,2*np.pi,M,endpoint=0)
dthetas =np.linspace(0,2*np.pi,1000,endpoint=1)

# values of the function at sample points
vals = fun(thetas)

A = np.zeros((M,2*N+1))

# construct the matrix of the system
A[:,0] = 1
for i in range(0,N):
    A[:,i+1] = np.cos((i+1)*thetas)
    A[:,N+i+1] = np.sin((i+1)*thetas)

coeffs = np.zeros(2*N+1)


B = (A.T)@A
print(np.shape(B))

# solve the system to find the coefficients
coeffs = np.linalg.solve(B,A.T@vals)

# function computing a trigonometric polynomial
def ptrig(x,c):
    res = np.zeros(np.shape(x))
    res = c[0]
    n = len(c)
    m = (n-1)//2
    for i in range(0,m):
        res = res+c[i+1]*np.cos((i+1)*x)+c[i+m+1]*np.sin((i+1)*x)
    return res

vals2 = ptrig(thetas,coeffs)

# plotting the result
plt.figure()
plt.plot(thetas,vals,'.b',label="Sample values")
plt.plot(dthetas,fun(dthetas),'g',label="Original function")
plt.plot(dthetas,ptrig(dthetas,coeffs),'r',label="Fitted trigonometric polynomial")
plt.legend()
plt.savefig("TrigPoly.png",dpi=100,bbox_inches='tight')
plt.show()

For the parameters chosen above the program outputs the following result. You can play with the input parameters to observe the changes.

Construct Pythagorean triangles with given constraints

This post will present some ideas related to the generation of all Pythagorean triangles satisfying a certain criterion, from an algorithmic point of view. Of course, there are infinitely many such triangles (integer sided with a right angle). Denoting by and the edges adjacent to the right angle and with the hypothenuse, we have the classical relation

When are all integers, it is possible to give a more precise result regarding the possible values of . Indeed, for any such triangle, there exist integers with , not both odd and coprime such that

These formulas are attributed to Euclid. The case corresponds to coprime and such a triangle is caled “primitive”. What is important to note is that this formula generates ALL pythagorean triangles exactly once, and this allows us to solve the following two questions:

1. Generate all pythagorean triangles with edges .

2. Generate all pythagorean triangles with or .

3. Generate all pythagorean triangles with .

Now let’s answer these questions one at a time:

1. In order to generate pythagorean triangles with edges we need to loop over all , coprime, not both odd with such that and for each primitive triangle, add to the list all its “multiples” by looking at the appropriate . One possible algorithm is

• set an empty list
• loop for to
  • loop for to , , not both odd
  • loop for to and add the triangle to the list
• in the end will contain the desired triangles

2. In this question one of the two legs or should be equal to . The only difficulty here is that we need to loop on the divisors of . Indeed, we have two cases:

• : for all divisors of , we should find the possible factorizations with not both odd and coprime, and then add the corresponding triangle to the list.
• : find all factorizations and check again that obtained are coprime and not both odd.

3. In this case . Therefore, one should loop on all divisors of and in each case solve

where , are coprime, not both odd. This can be done again with a loop.

These tools might come in handy when working on Project Euler problems, since often when dealing with integer sided quantities in a triangle, things can be reduced to pythagorean triangles. When you reach this step, it is enough to loop on these triangles and perform the requested operations.

Sum of floors of multiples of the Golden Ratio

Propose an algorithm for computing

for , where is the golden ratio. The sum is the sum of the first terms in the so called lower Wythoff sequence.

Solution: Computing a floor function and a multiplication is not complicated, therefore proposing a algorithm for computing is trivial. However, such an algorithm is not viable for .

The path to a much more efficient algorithm goes through the notion of Beatty sequence. For a positive irrational number the associated Beatty sequence is . This notion becomes interesting when noting the following fact. If is defined by then the Beatty sequences and cover all the positive integers. The latter sequence is called the associated complementary Beatty sequence. Two proofs of this fact can be found in the Wikipedia link above.

It is obvious now that is just the partial sum of the Beatty sequence associated to . Now let’s see what is the complementary sequence. A brief computation shows that the associated is given by , which shows that the terms which are not of the form are not far from it. In order to see how this can give an efficient algorithm, let’s follow the instructions below:

• denote by , the largest term in .
• looking at all the integers up to , in view of the result regarding the Beatty sequences, they are either of the form or of the form (in the associated complementary sequence).
• denote by the largest integer such that , which is given by . Then, it is not difficult to see that is the difference between the sum of all positive integers up to and .
• In the end we obtain

  

• by convention, state that or .

The last item above shows that a recursive algorithm can be implemented in order to compute . Moreover, since the computation of is reduced to the computation of where , the algorithm will converge very fast, since the sequence of upper bounds for converges exponentially to zero. Therefore, the algorithm obtained will have complexity and will work for extremely large values of , provided the language you are using can handle the expected precision.

The algorithm for computing shown above can be used, for example, in computing sums related to elements in the Wythoff table, which can be expressed using Fibonacci numbers and terms in .

Other values of , like for example (with ) lead to other types of sums for which the same type of algorithm can be applied. It is likely that even cases where is not explicitly related to through an integer may be handled through a similar procedure (to be confirmed).

Compute recurrent sequences fast

In many programming questions the computation of a recurrent sequence might be necessary. One of the famous such sequences is the Fibonacci sequence defined by and

Writing a code which computes for given is a no-brainer, but the complexity may vary much depending on the implementation. As it turns out, not everyone has the right reflexes when writing such a code. In the following I write some algorithmic variants and discuss their speed and complexity. In the end I’ll say a few words on how to compute quickly the -th sequence of any recurrence relation with fixed coefficients and fixed order.

First, remember that it is possible to have an exact formula:

where are solutions of and and are found from and . Note that are irrational, so when looking for an exact formula you might need to be careful for large . This could give you a algorithm. The formula is explicit, but you need to compute two exponentials so it’s not really free. Moreover, this shows that values of grow exponentially fast so you’ll quickly overflow if you’re not working in arbitrary precision. Moreover, if you need to compute for really large , this might not be the simplest approach, since you’re working with irrationals.

1. Recursion. The simplest way to code the computation of is using a recursive function. Create a function fibonacci which depending on the value of returns if or returns fibonacci()+fibonacci(). I remember using coding this when I was a beginner, thinking it was a nice idea… It turns out, however, that this approach has exponential complexity. If you’re not convinced, implement this and try to compute , etc. You’ll quickly see the exponential cost we’re talking about.

In order to see why, it is enough to note that when computing the algorithms needs to compute for all . Then for each one also needs to compute for all , etc. The only way this can be efficient if one stores the values already computed. This is called recursion with memoization, and should always be used when coding recursive tasks when there is a risk of computing the same value multiple times. However, the memoization has an memory cost.

2. Iterative sequence. The second approach is more efficient. You only need two variables, and which store and . Then at each iteration, store in a temporary variable, change into and store the value of in . It is obvious that this approach costs in execution time and in memory. If you need to store all values of , you will also have memory cost. As always, when doing this you can go up to in reasonable computing time (even higher, depending on your hardware, or programming language). However, reaching is generally out of the question.

3. Matrix exponentiation. One interesing approach, which is really fast is rewriting the second order recurrence as a first order matrix recurrence:

This implies that where . Therefore, in order to compute only a matrix exponential is needed. Using exponentiation by squaring this can be done in time. At the end you will obtain an (exact) integer result or a result modulo whatever integer you want. So this is the way to go if you need to compute for extremely large .

Please note that you can generalize this to recurrences of the form

It suffices to see that

so in the end you just need to compute the exponential of a bigger matrix. The complexity stays the same!

What to remember from this post?

• Don’t use recursion (without memoization) when computing recurrent sequences. You won’t get very far and you’ll waste your ressources!
• You can compute the -th value in recurrent sequences with constant coefficients with a complexity using matrix exponentiation. This works well even for huge when working modulo a fixed integer .

Are zero-order optimization methods any good?

The short answer is yes, but only when the derivative of gradient of the objective function is not available. To fix the ideas we refer to:

• optimization algorithm: as an iterative process of searching for approximations of a (local/global) minimum of a certain function
• zero-order algorithm: an optimization method which only uses function evaluations in order to decide on the next point in the iterative process.

Therefore, in view of the definitions above, zero-order algorithms want to approximate minimizers of a function using only function evaluations; no further information on derivatives is available. Classical examples are bracketing algorithms and genetic algorithms. The objective here is not to go into detail in any of these algortihms, but to underline one basic limitation which must be taken into account whenever considering these methods.

In a zero-order optimization algorithm any decision regarding the choice of the next iterate can be made only by comparing the values of  at different evaluation points. For example, look at the classical trisection algorithm for minimizing a unimodal function, i.e. a real function defined on which is decreasing on and increasing on :

• given the bracketing choose the points and .
• compare the values and in order to decide on the next bracketing interval:if then

  if then

• stop when the difference is small enough.

One question immediately rises: can such an algorithm reach any desired precision, for example, can it reach the machine precision, i.e. the precision to which computations are done in the software used? To fix the ideas, we’ll suppose that we are in the familiar world of numbers written in double precision, where the machine precision is something like . I will not go into further details regarding this, since people do whole courses on this. Note that Matlab, Octave and Numpy are well known languages in which the default setup uses this machine precision.

More precisely, the real numbers are stored as floating point numbers and only digits are relevant in the computations. When adding two numbers whose ratio is smaller than the result will be . This is due to the fact that adding the numbers would require us to shift the decimal point to the same position. However, since the ratio is smaller than , when shifting the decimal point to the same position, the resulting number will contain zeros on the first or more significant digits in . Therefore the addition will not change any significant digits in .

This issue related to computations done in floating point precision makes the question of comparing and in the algorithm above pointless when and are close to the minimum. In order to identify the source of the potential problem let us write the Taylor expansion of the function around the minimum . First, note that at the minimum, which leaves us with

where denotes higher order terms. You may note that neglecting higher order terms shows that the function looks like a quadratic function near the optimum. This is a good thing to keep in mind and students should not hold any grudge for professors who illustrate the behavior of optimization algorithms on quadratic functions. As it turns out, any good algorithm needs to behave well for quadratic functions, since every function, near a minimum is eventually quadratic.

Coming back to our previous paragraph, it turns out that the computer won’t be able to tell the difference between and if

where is the machine precision. Therefore, if

the machine won’t tell the difference between and . Since the second factor in the above expression is well defined and positive in most cases, it turns out that we will not be able to decide if we are closer than to the solution based only on the values of the function .

Therefore, when using a zero-order algorithms for a function with you cannot decide based only on the function values if you are closer to the optimum than the square root of the machine precision . This is quite a limitation, and in practice it should save lots of pointless function evaluations.

When using derivatives or gradients there is no such problem, since we can decide if we are close enough to by comparing the derivative to zero, and this comparison is valid up to the precision to which we can compute .

In conclusion, use zero-order methods only if gradient based methods are out of reach. Zero-order methods are not-only slowly convergent, but may also be unable to achieve the expected precision in the end.

On Heron’s algorithm – approximating the square root

Heron’s algorithm computes an approximation of the square root of a positive real number using the iteration

1. Show that the recurrence defined above is Newton’s algorithm applied for finding a zero of the function . Show that the sequence converges to when is close enough to . What is the order of convergence of the sequence?In the following we consider the following error measure:

   

   This choice simplifies the computations and is an estimation of half the relative error when .

2. Show that the sequence verifies

   

   Deduce an explicit formula for in terms of . Show that the sequence converges to for every choice of the initial condition .

3. Show that for every there exists an integer such that . Write a modified Heron algorithm which reduces the computation of the square root of to the computation of the square root of a real number .The choice of the initialization determines and has an influence on the speed of convergence. In the following we restrict ourselves to the case and we look at the choice of the initial condition. We will suppose that the initial condition depends on and we will consider the following cases: , a constant function and , a polynomial of degree in . In every case we will try to minimize for .
4. Show that if then is minimal for . Find the explicit value of . (Hint: You could represent graphically the function .)
5. Suppose that is a polynomial of degree in . We want to find the coefficients and which minimize with .
   (a) Prove that the maximum of is found for .
   (b) Show that in order to minimize it is necessary that and that .
   (c) Conclude that the and which minimize are . Find an exact expression for .
6. For compare the number of iterations necessary to arrive at a demanded precision ( for the two initializations studied: and with . Comment on the interest of using the initialization of degree compared to the constant initialization
7. Evaluate the number of operations used for the two cases and conclude which of the two approaches is more advantageous.
8. Try to find the initialization of second degree which is optimal and perform a similar investigation to see if from the point of view of the number of operations this choice is well adapted or not.

p-Laplace equation – Fixed point approach

The -Laplace problem is the generalization of the usual Laplace problem and is defined for in the following way

We are not going into details regarding the functional spaces corresponding to and since in the end we are interested in showing a way to solve this equation numerically.

Project Euler tips

A few years back I started working on Project Euler problems mainly because it was fun from a mathematical point of view, but also to improve my programming skills. After solving about 120 problems I seem to have hit a wall, because the numbers involved in some of the problems were just too big for my simple brute-force algorithms.

Recently, I decided to try and see if I can do some more of these problems. I cannot say that I’ve acquired some new techniques between 2012-2016 concerning the mathematics involved in these problems. My research topics are usually quite different and my day to day programming routines are more about constructing new stuff which works fast enough than optimizing actual code. Nevertheless, I have some experience coding in Matlab, and I realized that nested loops are to be avoided. Vectorizing the code can speed up things 100 fold.

So the point of Project Euler tasks is making things go well for large numbers. Normally all problems are tested and should run within a minute on a regular machine. This brings us to choosing the right algorithms, the right simplifications and finding the complexity of the algorithms involved.

ICPC 2015 World Final Problem A

This is the solution to Problem A from the International Collegiate Programming Contest. The list of problems can be found here. This first problem consists simply of reading the parameters of the function defined below, and computing its values on the set . Then, you need to find the maximum decrease in the function values.

The inputs are parameters , the function is

and the values to be considered are . Below is a Matlab code which works, at least for the given Sample Inputs. This should be optimized by removing the for loop. It is instant for but it should be modified to work until .

function res = Prob1_2015(p,a,b,c,d,n)
dis = 1:n;

vals = p*(sin(a*dis+b)+cos(c*dis+d)+2);

mat = zeros(n,1);
for i = 1:n
mat(i) = vals(i)-min(vals(i:end));
end

res = max(mat)

The sample outputs are:

 >> Prob1_2015(42,1,23,4,8,10)  
      104.855110477394

 >> Prob1_2015(100,7,615,998,801,3)
      0

 >> Prob1_2015(100,432,406,867,60,1000)
      399.303812592112

New version which works for large . It suffices to eliminate the computation of the min at each of the phases of the iteration. This solves the problem in 0.2 seconds for

function res = Prob1_2015(p,a,b,c,d,n)
dis = 1:n;
vals = p*(sin(a*dis+b)+cos(c*dis+d)+2);
mat       = zeros(n,1);
lastmax   = vals(1);
for i = 1:n
  lastmax   = max(lastmax,vals(i));
  mat(i)    = lastmax-vals(i);
end
res = max(mat)

Identifying edges and boundary points – 2D Mesh – Matlab

A triangulation algorithm often gives as output a list of points, and a list of triangle. Each triangle is represented by the indexes of the points which form it. Sometimes we need extra information on the structure of the triangulation, like the list of edges, or the list of boundary points. I will present below two fast algorithms for doing this.

Finding the list of edges is not hard. The idea is to go through each triangle, and extract all edges. The algorithm proposed below creates the adjacency matrix in the following way: for each triangle we set the elements (and their symmetric correspondents) to be equal to .

In order to find the points on the boundary (in two dimensions), it is enough to look for the edges which are sides to only one triangle. We can do this using the adjacency matrix. Note that if is the adjacency matrix, then stores the number of paths of length (two sides) between two points of the triangulation. Note that any edge which is not on the boundary will contain the starting and ending point of two paths of length . If is a triangle such that points are on the boundary, then (there is one path of length going through . We also have . Conversely, if and then are connected, and there is a unique path of length going from to . Thus, is an edge on the boundary. Therefore, we just need to identify the indexes such that there exists with .

Below are two short Matlab codes doing these two algorithms. I guess they are close to being optimal, since only sparse and vectorized operations are used.

%p is the list of points
%T is the list of triangles, ap is the number of points
%this computes the adjacency matrix
A = min(sparse(T(:,1),T(:,2),1,ap,ap)+sparse(T(:,2),T(:,3),1,ap,ap)+sparse(T(:,3),T(:,1),1,ap,ap),1);
A = min(A+A',1);
% this finds the boundary points, whose indexes are stored in Ibord
B = A^2.*A==1;
Ibord = find(sum(B,2)>0);

Distance from a point to an ellipsoid

Here I will discuss the problem of finding the point which realizes the minimal distance from a given to an ellipsoid. First I will look at the theoretical properties which will lead us to and then I will present a numerical algorithm to find the nearest point.

This problem is interesting in itself, but has some other applications. Sometimes in optimization problems we need to impose a quadratic constraint, so after a gradient descent we always project the result on the constraint which is an ellipsoid.

Let’s begin with the problem, which can be formulated in a general way like this: we are given , and a matrix , positive semidefinite, which gives an ellipsoid . Without loss of generality assume that is diagonal (if not, we diagonalize it using an orthogonal matrix and work in a new coordinate system). The problem becomes

(note that the square is not important here, but simplifies the computations later)