Analyzing lfo-run

I have written tests that compare lfo-run from TISEAN to the ported version lfo_run. The test that uses amplitude.dat works perfectly, but when I analyzed the results both programs/function gave for henon (Henon Maps) I ran into some problems. I will attempt to describe them.

Input data

The problems occur when analyzing a 1000 element Henon map (henon(1000)). For all of the implementations if I used a simple call with default parameters (m = 2, d= 1) the programs would quit due to a matrix singularity. The problems arose when (m = 4, d= 6) was used. With these parameters the program gave various results for various implementation methods.

It is important to note that the prediction that I was testing tried to predict 1000 future elements (default for all implementations) on the basis of given 1000 elements.

Implementations

There are 3 implementations I used:

1. The TISEAN implementation (uses lfo-run)
2. The implementation similar to 1. but compiled as c++ and wrapped in enough code to run as m-file (uses __lfo_run__ and invert_matrix())
3. The implementation that uses Matrix::solve() method

I tried to find out if maybe method 1. and 2. do not differ due to a bug that was introduced while porting. I therefore ported it twice (the second time to a very rudimentary stage) and both times the same results were encountered. I do not understand why there is a discrepancy between between these two implementations.

Discrepencies

Since the goal of this project is to port TISEAN functions I have compare implementation 1. with 2. and 1. with 3, to see what differences I come across.

A comparison between implementation 1. and 2. results in an error from implementation 2. The function generates about 700 elements (of the default 1000) and then throws an error that the forecast has failed.

The comparison between implementation 1. and 3. is much more fruitful as the results are the same for about 150 elements and then they begin to differ (see Fig. 1.)

Fig. 1 Comparison between the TISEAN implementation and using Matrix::solve() in TISEAN package from Octave

These results can be achieved by cloning the repo, doing make run and running the script:

cd tests/lfo_run/; test_lfo_run



Analysis

Before I give my suggestions for what is the cause for these discrepancies I would like to discuss another interesting discrepancy. This discrepancy is the maximum difference between the solution of the equation system obtained from implementation 2. and 3. When using implementation 3. for the forecast this value was 8.5e-14, but when using implementation 2. for the forecast this difference was 5e-13.

I believe this is because the computational error is accumulated throughout the program. Each new forecast point is dependent on the previous ones. Moreover the Kahan algorithm (compensated summation) is never used in the TISEAN implementations. Even matrix multiplication (as seen e.g. in multiply_matrix())  uses the simple, but error accumulating for(...) sum += vec[i].

As to why implementation 1. and 2. give different results I have two theories: either there still is a bug which I was unable to detect, or some compilation difference (e.g. a linked library) between the TISEAN program written in C (lfo-run) and the TISEAN package function written in C++(__lfo_run__).

Summary

The question that poses itself is whether this warrants a rewriting of other TISEAN function that use the simple summation, or if this problem can be ignored. The authors of TISEAN said in their introduction that blindly using the programs they wrote may result in unintended or even wrong results. Trying to predict 1000 elements of a 1000 element Henon Map using first order local linear prediction might be considered a bad use-case.

Progress report

Since the last post I wrote a tutorial on http://wiki.octave.org/TISEAN_package. I also ported:

• lzo_gm
• lzo_run
• ikeda
• lfo_run
• lfo_ar
• lfo_test
• rbf
• polynom

During my work I discovered that polynom has similar functions (polynomp, polypar, polyback) which provide extra options for performing polynomial fits. I will not include those functions in the project now, but they have a high priority once I finish all of the functions I outlined for this project.

These newly ported functions aren't completely polished (some need demos and documentation) but they pass tests and don't have memory leaks. Once I clean these functions up and port xzero the last function in this section, I intend to create version 0.1.0 of the package. With this version I intend to branch out the repo to have a 'devel' and a 'stable' branch.

Afterwards, I will add more information to the tutorials on the wiki page.