Examination projects

The algorithm to pick a project:

• pick the project with your ordinal number from the spreadsheet (check it once again);
• if your number, n, is larger than the number of projects, N, continue counting from the first project (that is, if n>N, your exercise is the remainer of integer division of n by N);

You have the right to choose instead one of the more complicated --- in my opinion --- voluntary project below. The above algorithm for picking your project applies also to voluntary projects.

Remember to create a readme.txt file with a description of what you were doing.

time now: 01/Apr/2026-20:19
next run: 02/Jul/2013-10:00

Ok, it is reshuffled for 02/Jul/2013-10:00 run: go ahead and good luck... :)

1. Multi-dimensional pseudo-random adaptive recursive (non-stratified: simply subdivide in all dimensions) integrator in a rectangular volume.

   Do not waste points.

   Check how far in the number of dimensions you can go without overwhelming the computer (most probably running out of stack, but in any case up to 100M allocated memory and up to 5min nice run time).

2. Two-dimensional integrator for integrals in the form
   ∫abdx ∫d(x)u(x)dy f(x,y)
   which consecutively applies your favourite adaptive one-dimensional integrator along each of the two dimensions. The signature might be something like

   double integ2D(double a, double b, double d(double x), double u(double x), double f(double x, double y), double acc, double eps)
   

3. Quantum annealing (minimization).

4. Inverse iteration algorithm for eigenvalues (and eigenvectors).

5. One-dimensional adaptive integrator which at each iteration subdivides the interval not into two, but into three sub-intervals.

   Use a closed classical quadrature unless the integrand returns INFINITY in which case switch to an open quadrature.

   Compare with your other integrators.

6. Multidimensional pseudo-random integrator -- that is, plain Monte Carlo: Q(N)=V⋅average({f(x_i)}), ΔQ(N)=V/√N⋅variance({f(x_i)}) --- with asymptotic error estimate:
   • For a set of several different numbers of sample points, N₁,..., N_n (for example, 1000, 2000, 3000) estimate --- using plain Monte-Carlo --- the integral, Q_i=Q(N_i), and the error (through the variance, as usual), ΔQ_i=ΔQ(N_i). Of course you should do it in one run, storing the data-points (N_i, Q_i, ΔQ_i) along the sampling process.
   • Fit this data set with error-bars, {(N_i, Q_i, ΔQ_i) | i=1,...,n}, using the fitting function f(N)=c₀+c₁/√N. The coefficient c₀ and the corresponding error Δc₀ are now the asymptotic estimates of the integral and the error.
   • If the error is yet too big, add an extra data-point {(N_n+1, Q_n+1, ΔQ_n+1)} to your data and redo the fit; repeat until the error is small enough.

   Check whether this algorithm gives any improvement over plain Monte Carlo.

   Try use the fitting function f(N)=c₀+c₁/√N+c₂/N.

7. Multidimensional quasi-random integrator with asymptotic error estimate:

   • Throw N quasi-random points and estimate, say, Q(1/2 N), Q(2/3 N), Q(5/6 N), Q(N) (where Q(n) is the quasi-random estimate of the integral using n points).

   • Make an ordinary least-squares fit to these data-points using the fitting function f(N)=c₀+c₁/N.

     (Although the integration errors with quasi-random numbers are not readily available, you can still promote the contributions from the estimates with larger samples by introducing artifitial errors inversly proportional to the number of sampling points.)

     The coefficient c₀ is apparently the asymptotic estimate of the integral (under the assumption that the integration error behaves as O(1/N)).

   • Add some more points and repeat the asymptotic estimate. The difference between the old and the new c₀ gives an estimate of the error.

   • Repeat until the given accuracy goal is achieved.

   Try fitting more data-points; try using more term in the fitting function, say f(N)=c₀+c₁/N+c₂/N².

8. Quasi-Newton's minimization with SR1 update.

9. Implement a (sub-)spline of a data set {x_i, y_i| i=1,..,n; p₁, p_n} --- where the user also supplies the derivatives p₁ and p_n at the end-points x₁, x_n --- using the following algorithm:
   • for each inner point x_i, i=2,..,n-1, build a cubic interpolating polynomial through the points x_i-1, x_i, x_i+1 and evaluate, using this polynomial, the first derivative, p_i, at the point x_i;
   • now you have a data set {x_i, y_i, p_i| i=1,..,n} where the function is tabulated together with its derivative: build a cubic sub-spline for it,
     S_i(x)= y_i +b_i(x-x_i) +c_i(x-x_i)² +d_i(x-x_i)³,
     where for each interval the three coefficients b_i, c_i, d_i are determined by the three conditions, S_i(x_i+1)=y_i+1, S'_i(x_i)=p_i, S'_i(x_i+1)=p_i+1. (See the subsection 'Akima sub-spline interpolation' for the inspiration)

   Check whether it makes smaller wiggles than quadratic and cubic splines.

   Now, if the user supplies NAN (from math.h) for one or both of the end-point derivatives, you have to estimate p₁ and/or p_n from the the same cubic interpolating polynomial that you used to determine correspondingly p₂ and/or p_n-1.

10. Symmetric row/column update of a size-n symmetric eigenvalue problem.

    The matrix to diagonalize is given in the form

    A = D + e(p) u^T + u e(p)^T

    where D is a diagonal matrix with diagonal elements {d_k, k=1,...,n}, u is a given column-vector, and the vector e(p) with components e(p)_i=δ_ip is a unit vector in the direction p where 1≤p≤n.

    Given the diagonal elements of the matrix D, the vector u, and the integer p, calculate the eigenvalues of the matrix A using O(n²) operations.

    Although this problem can be solved through two rank-1 updates (see the notes), you probably want to use the following method which should be simpler.

    Without loss of generality we can assume that the component u_p is equal zero (if it is not, we can make it so by simply redefining the diagonal element d_p→d_p+2u_p, u_p→0).

    The eigenvalue equation for the updated matrix reads

    (D - λ) x + e(p) u^T x +u e(p)^T x = 0,

    where x is an eigenvector and λ is the corresponding eigenvalue.

    The component number p of this equation reads

    (d_p - λ) x_p + u^T x = 0,

    while the component number k≠p reads

    (d_k - λ) x_k + u_k x_p = 0.

    Dividing the last equation by (d_k-λ), multiplying from the left with ∑_ku_k, substituting u^T x from the first equation and dividing by x_p gives the sought secular equation,

    -(d_p-λ) + ∑_k u_k²/(d_k-λ) = 0.

    The roots of this equation give the updated eigenvalues (check the derivations thoroughly though, it is only an idea, I didn't try it myself).

    This equation is similar to the secular equation of rank-1 update (see the notes) and has similar properties.

11. Symmetric rank-1 update of a size-n symmetric eigenvalue problem.

    The matrix A to diagonalize is given in the form

    A = D +uu^T,

    where D is a diagonal matrix and u is a column-vector.

    So, given the diagonal elements of the matrix D and the elements of the vector u find the eigenvalues of the matrix A using only O(n²) operations (see section ''Eigenvalues of updated matrix'' in the book).

12. Implement a subroutine that calculates natural logarithm of complex argument z via the integral

    ln(z) = ∫1z dt/t

    by directly calculating the integral of the complex-valued integrand of the complex argument along a straight line in the compex plan.

    Think, what happens if you integrate not along a straight line but rather make a circle (or two) around zero :)

13. Multidimensional pseudo-random (plain Monte Carlo) vs quasi-random (Halton and/or lattice sequence) integrators: investigate the convergence rates (of some interesting integrals in different dimensions) as function of the number of sample points.

14. Implement a subroutine that calculates the error function of complex argument z,

    erf(z) = 2/√π ∫0z exp(-t2)dt ,

    by directly calculating the integral of the complex-valued integrand of the complex argument along a straight line in the compex plan. You can use complex type (C, C++); complex<double> (C++); and python's module 'cmath'. Try different recursive adaptive integrators and pick the most effective one for the job. Investigate the realistic accuracy of this method of calculating the error function; the realistic region in the complex plane where you can reliably calculate the error function. Compare with the built-in error function from math.h.

15. Akima sub-spline.

16. Build a cubic (sub-)spline, S_i(x)= y_i +b_i(x-x_i) +c_i(x-x_i)² +d_i(x-x_i)³, for a data set, {x_i, y_i, y'_i| i=1,..,n}, where the function is tabulated together with its derivative (see the subsection 'Akima sub-spline interpolation' for the inspiration) -- for example, as a result of a numerical integration of an ordinary second order differential equation. Now make the second derivative of the spline continuous by increasing the order of the spline to four, for example, in the form S_i(x)= y_i +b_i(x-x_i) +c_i(x-x_i)² +d_i(x-x_i)³ +e_i(x-x_i)²(x-x_i+1)² .

17. One-dimensional adaptive integrator which initially uses a closed embedded classical quadrature (allegedly more precise than an equivalen open), but switches to open quadratures on intervals where it encounters INFINITY from the integrand.

    Improvement: on intervals with INFINITY switch to Clenshaw-Curtis variable substitution algorithm.

18. Quasi-Newton's root finding with Broyden's update

19. ODE driver that checks that |δyk| < (δ+ε|yk|)√(h/(b-a)) where k=1...n counts the components of the vector-function y (that is, instead of the simpler condition ||δy|| < (δ+ε||y||)√(h/(b-a))) together with a two-and-one-third-step stepper (which uses one previous point and makes and auxiliary h/3 step).

20. QR-decomposition using Givens rotations.

21. Bi-linear interpolation on a regular grid in two dimensions. The signature of the interpolating subroutine can be

    double bilinear(int nx, int ny, double* x, double* y, double** z)

    or

    double bilinear(gsl_vector* x, gsl_vector* y, gsl_matrix* z)

22. Quasi-Newton's minimization with Broyden's update

Voluntary projects

1. Evolutionary minimization algorithm of your choice.

2. Few smallest eigenvalues of a matrix using Lanczos tridiagonalization plus recursive diagonalization (amounts to rank-1 update).

3. Multi-dimensional Monte Carlo integration with recursive stratified sampling. (Try use the largest |⟨f⟩left-⟨f⟩right| as the criterion for choosing the bisection direction instead of the largest max(varianceleft,varianceright).)

4. Simulated annealing with downhill simplex or your other favourite local minimizer.

The functions you could use to test your 1D integrators

• ∫01 exp(x) dx
• ∫01 ln(x) dx
• ∫01 x1/2 dx
• ∫01 x3/2 dx
• ∫01 x-1/2 dx
• ∫01 (x>0.3?1:0) dx
• ∫03 int(exp(x)) dx