The Pain of Java Matrix Libraries

Looking for a good Java Matrix (and actually also math) library, I was a bit surprised to find out there does not seem to be any really serious one still maintained.

Sure, there is Apache Commons Math, but it is still changing a lot, and it is not very performance optimized yet, while it has been active for several years already. There is also Java3D, it does Matrix through GMatrix, but not much linear algebra and if you look at their implementation, it is very basic, not performance oriented.

The other candidates seem to be 1-man projects that can disappear any other day (some of them look quite good like ojalgo, most of them are not interesting). Then you also have the serious but not maintained anymore Cern Colt library.

Compared to C/C++, the Java world is worrying if you want to do maths.

In those libraries, a dense matrix of double can be implemented two ways:

  • by maintaining internally a double[][]. Usually those libraries allow for not copying the array, so it can be neat if your interfaces have this kind of arrays.
  • by maintaining internally a double[]. The reason is for performance, but then each time you build a matrix from a double[][], an expensive copy will happen. So you need to use the Matrix object in your interfaces instead of double[][].

This is a pain because you can be very quickly stuck in one or the other Matrix library. A “solution” is to have your own interface, but that is a pain to write. There is UJMP, but it can hide some important methods (like transpose and multiply in one go from Colt or the ability to reuse an existing matrix in various operations to avoid allocating new memory), it is a students project (like parallel colt), but if it was a standard, it could be much more interesting.

In summary it does really look like scientific people, universities don’t use Java for computation otherwise Colt surely would have been maintained.

Java Calendar Is Broken On JVM Upgrade

We ran into an interesting issue with TimeZone and Dates. If you print the same date on different JVMs, it might show a different printed date.

The reason behind this is the daylight saving time conventions. An old JVM won't necessarily have the same daylight saving time for a given TimeZone than a latest JVM, and therefore will interpret the date differently.

Here is the output of a very simple program on 2 different JVMs. The number is the long number used to represent the date. I only use SimpleDateFormat with different TimeZone:

JVM 1.5.0_12
loading date=Sat Jul 18 06:59:36 CEST 2009, number=1247893176505
using formatter in EST: 7/17/09 11:59 PM
using formatter in Asia/Singapore: 7/18/09 12:59 PM

JVM 1.5.0_20
loading date=Sat Jul 18 06:59:36 CEST 2009, number=1247893176505
using formatter in EST: 7/18/09 12:59 AM
using formatter in Asia/Singapore: 7/18/09 12:59 PM

The source code:

public class DateBug {

private static String FILE_NAME = "datebug.txt";

public static void load() throws IOException {
FileReader fr = new FileReader(FILE_NAME);
BufferedReader br = new BufferedReader(fr);
String l = br.readLine();
br.close();
long time = new Long(l);
Date d = new Date(time);
System.out.println("loading date="+d+", number="+d.getTime());
SimpleDateFormat formatter = new SimpleDateFormat();
formatter.setTimeZone(TimeZone.getTimeZone("EST"));
System.out.println("using formatter in EST: "+formatter.format(d));
formatter.setTimeZone(TimeZone.getTimeZone("Asia/Singapore"));
System.out.println("using formatter in Asia/Singapore: "+formatter.format(d));
}

public static void saveNew() throws IOException {
Calendar c = Calendar.getInstance(TimeZone.getTimeZone("EST"));
c.set(2009, 06, 17, 23, 59);
Date d = c.getTime();
System.out.println("saving date="+d+", number="+d.getTime());
SimpleDateFormat formatter = new SimpleDateFormat();
formatter.setTimeZone(TimeZone.getTimeZone("EST"));
System.out.println("using formatter in EST: "+formatter.format(d));
formatter.setTimeZone(TimeZone.getTimeZone("Asia/Singapore"));
System.out.println("using formatter in Asia/Singapore: "+formatter.format(d));

FileWriter fw = new FileWriter(FILE_NAME);
PrintWriter pw = new PrintWriter(fw);
pw.println(d.getTime());
pw.close();
}

public static void main(String[] args) throws IOException {
System.out.println("JVM "+System.getProperty("java.version"));
if (args.length == 1) {
if (args[0].equals("save")) {
saveNew();
}
} else {
load();
}
}
What does this mean? This means that if you entered in a GUI in the first JVM a particular date & time using EST time zone. This will change when you read back in the second JVM.
This suggests that if you want to keep the same dates, you are better off saving in UTC where daylight saving time is not used and trick DateFormat. But I have to say this looks quite ugly.

Implicit Finite Differences Method For Pricing Barrier Option


While trying to price a simple knock down and out barrier option, I encountered several difficulties I did not expect with the implicit finite differences method. The explicit method has less issues with barrier options pricing. I will show here what the tricky parts are and why explicit seems simpler in this case.

The full article is here (pdf) or here (html) (the later is not very well formatted).

Algorithms used in this article can be found at http://code.google.com/p/javamc/

Pulseaudio Nightmares - Pure ALSA to the Rescue

In the latest stable Ubuntu (9.04), pulseaudio still does not work reliably on my hardware (intel HDA with digital SPDIF out). I upgraded to 9.10 and had even more problems. Sound always worked on boot, but often broke down after a while. And I could not find easy ways to make it work other than rebooting... Killing/restarting pulseaudio, looking at the processes using snd did not work.

One thing works wonderfully, pure ALSA. To have multilple apps sharing ALSA, I just use dmix. As I use digital out, there is no mixer, but ALSA can provide one through softvol. It works really well. ALSA is already not that simple to configure/setup properly, but with pulseaudio on top, welcome to your worst configuration nightmares.

Here is the .asoundrc I use:
 
pcm.amix {
type dmix
ipc_key 50557
slave {
pcm "hw:0,1"
period_time 0
period_size 1024
buffer_size 8192
}
bindings {
0 0
1 1
}
}

pcm.softvol {
type softvol
slave {
pcm "amix" #redirect the output to dmix (instead of "hw:0,0")
}
control {
name "PCM" #override the PCM slider to set the softvol volume level globally
card 0
}
}


pcm.!default {
type plug
slave.pcm "softvol" #make use of softvol
}

Java int Overflow Behavior

A coworker recently asked me if there was a guaranteed behavior in the case of int overflow. He gave the specific example on:

can we rely that int x = Integer.MAX_VALUE + 1 is the same for every JVM on any platform?

I thought the answer would be easy to find in the Java specifications document. But I was wrong. It is not clearly defined.

I found a trick that suggests this behavior is indeed standard and will stay the same, it is related to type casting. Java guarantees that the cast of a long just truncates the long to int precision. Therefore if

long l = Integer.MAX_VALUE;
l = l +1;
x = (int)l;

x has a guaranteed value. http://java.sun.com/docs/books/jls/third_edition/html/conversions.html paragraph 5.1.3.

Static Fields and Inheritance

Someone asked me recently to find out the real reason why the code from this thread fails. This is a fairly bad code, and not even a very good way to point out the problem. But the question is nonetheless interesting.

class Toto extends TotoParent{

final static Toto a = new Toto ("a");

public Toto(String a){
super(a);
}
}

import java.util.ArrayList;
import java.util.List;

public abstract class TotoParent {

static List list = new ArrayList();

public TotoParent(String a) {
list.add(a);
}

protected static List get() {
return list;

}
}

import org.junit.Test;
import static org.junit.Assert.*;

public class TotoTest {

@Test
public void testGet(){
assertEquals(1, Toto.get().size());
}
}
I am quite used to static initialization, and would have answered the same as the first answer in the thread:
"Get is static and associated with TotoParent, so that is the same as calling TotoParent.get().size()". I would have even thought that the compiler would compile the call Toto.get() to TotoParent.get(). But running javap, you can see it is still compiled as TotoParent.get(). So there is still a lookup done. This is why the first answer is actually not that correct.

The important bit here is that Toto is never initialized, even if we call Toto.get(). The java specs (invaluable reference) explains clearly that calling a static method not declared in the class does not initialize the class.

Calling Toto.get() is not exactly the same as calling TotoParent.get().
If TotoParent.get() called another TotoSuperParent.get():
Toto.get() -> TotoParent.get() -> TotoSuperParent.get()
We compile then later we change to make TotoParent have a specific implementation of get(). Toto will then be automatically aware of it, without even recompiling it.

http://java.sun.com/docs/books/jls/third_edition/html/execution.html
paragraph 12.4.1

Benchmarking Languages Is Difficult

I often looked at the famous computer languages shootout for fun. Recently I noticed they had the infamous thread ring test. I posted not very long ago several blog entries about it showing how silly this test was.

Looking at the existing Java implementation for the test I decided to try to submit the tricky one using a pool of thread, and pooling message processing rather creating 1 thread per node. To my surprise, it was accepted without questions and I did have the best score for a Java program for a while. Shortly after someone else copied my program and got rid of various stuff not useful for the particular benchmark (breaking the interesting part of the design) and got accepted as well with of course a better result.

I decided to see if I could make an even more silly program - tailored for the test only. I managed to be orders of magnitude faster - 1 thread, no synchronization, everything processed in a FIFO (linkedlist) queue. This is actually a standard way to reimplement recursion. But I was honest enough not to hide that I consider that kind of program to cheat the test and got my entry in the "interesting alternatives".

In reality there is no difference in the "cheating" between my new program and the program that got accepted in the official list, they both cheat by using only 1 thread and process everything 1 by 1. There is not 1 thread per node in any of the program, and they can avoid any concurrency issues. One "looks" better because it uses a pool of 503 threads (but really use only 1 or 2 threads) and the other does not hide its use of 1 thread for processing. But this is not evident to people accepting the programs.

When I look at the haskell code, I can not really tell if it is creating 503 threads in the language or a pool or ..., you have to know each language quite well and sometimes it is not that easy to define what cheating is. Therefore this kind of benchmark is a bit disappointing. One should force the use of the same algorithm. But can you do so (a functional language won't use the same algo as a procedural one)?

Cholesky & Jakarta Commons Math

In Finance, Cholesky is a useful way to decompose Matrix. It is not so simple to find a BSD licensed code using cholesky (most of them are GPL like this one). There is one in Apache Commons Maths library, which is a very interesting library. However for performance, it is still not very practical for some things like Cholesky.

Looking at the source one can easily understand why. I did a small (many people will say not representative 1 million loop test) and finds out:

  • cholesky GPL= 5.4ms
  • cholesky BSD=37.1ms

So BSD code is 7 times slower! Of course it can do a bit more and has many checks of validity, but still. It shows it is not easy to do Math libraries, because some people will care a lot about this performance difference, and some other people won’t but will like the other “features”.

Hull American Option Price Fallacies


Hull says American put is best exercised immediately and american call is optimal at expiry like a european. Is this really true?

At first it seems really clever and model show clearly this. But if we change the market assumptions only a tiny bit, everything falls down.

I could not detail everything in a blog post so I created a static web page about it. Everything was produced in Java using algorithm found in popular books and graphs through JFreeChart.

On Quasi Random Numbers - MersenneTwister vs Sobol precision in Asian Option Pricing

While starting a side project that does Monte Carlo pricing in Java (http://code.google.com/p/javamc/ - nothing yet there I am waiting for Mercurial repository support), I wondered what was the importance of quasi random numbers versus more regular pseudo random numbers in Monte Carlo simulations.

This brought me to read more carefully several books about Monte Carlo and Finance (Haug Option Pricing, Sobol Primer on Monte Carlo, and Glasserman Monte Carlo Methods in Finance Engineering). I had quite a hard time to understand why the dimension of the quasi random generator was so important to price an asian option. Intuitively I thought the averaging points of an asian option were all on the same path, so they should be using the same random generator. This is very wrong as one does not care about the path in the first place but just in simulating each point in the average (using the regular black and scholes hypothesis). Finding the estimation for the average on the given points forces to use independent random generators for each point, because we want to approximate the estimation by the sum over those random points for each point.

There is another simple argument to explain why independence of the random generators is so important. If we use the same generator for each point, then each point will move exactly the same way at each simulation. The average of those point will therefore behave exactly the same way as if there was only 1 point using the same generator. And we don't price an asian anymore but just a regular vanilla option.

Using a pseudo random generator, one does not see the problem of dimension, because we can create N independent dimensions by just taking numbers N by N on a pseudo random generator. So effectively having 1 or N dimensions is the same on a pseudo random generator.

Still I wrote a small test to see if a 1D quasi random generator was so bad when simulating N dimensions (taking values N by N on the quasi random generator). Here are the results:

MersenneTwister vs MersenneTwister on 10D asian:
14:43:51,111 INFO MonteCarloSimulationTest:114 - 867970000 -- expPrice=0.978958644504466
14:43:51,428 INFO MonteCarloSimulationTest:120 - 314619000 -- expPrice=0.9733220318545934
14:43:51,430 INFO MonteCarloSimulationTest:122 - relative difference=-0.005757763804951897
can be as high as 2%

Sobol vs MersenneTwister on 10D asian:
14:48:46,909 INFO MonteCarloSimulationTest:115 - 980209000 -- expPrice=0.9895032774079221
14:48:47,345 INFO MonteCarloSimulationTest:121 - 433685000 -- expPrice=0.9790264042895171
14:48:47,348 INFO MonteCarloSimulationTest:123 - relative difference=-0.010588012548932534
about 1% it is actually bounded by MersenneTwister precision.

Sobol vs Sobol1D on 10D asian:
14:47:08,614 INFO MonteCarloSimulationTest:115 - 717444000 -- expPrice=0.8810736428068913
14:47:08,925 INFO MonteCarloSimulationTest:121 - 308499000 -- expPrice=0.9791449305055208
14:47:08,927 INFO MonteCarloSimulationTest:123 - relative difference=0.11130884290920073
about 10% and stays that way even when increasing number of simulations.


Using an asian rate with 10 points, we see that Sobol1D will always give a very bad estimate, no matter the number of simulations. While Sobol used properly will give (much) better precision for less iterations. So even though there is the word random in quasi random, the numbers are very far from being random or even behaving like random numbers. It helped me to read about Van der Corput and Halton numbers to really understand quasi random numbers.

Previous

Next