Zachary Danziger

Lingji,

I cannot see a way to calculate the coupling distance during the recursive calls to c(i,j) without additional overhead. I feel that adding computational cost is not justified because, in general, the coupling sequence is not unique and therefore is not especially informative. The coupling sequence is just any allowable (i.e., follows the forward movement rules) sequence between points on P and Q that never has distance greater than cm. In the typical case there are a great many allowable sequences because interchanging many points with small distances does not affect cm if there is a much larger distance later in the sequence.

Here is my compromise: if the user requests the coupling sequence then I calculate the cm as usual, and at the end of the code I loop through the CA variable choosing one workable coupling sequence. I hope this pleases the fans.

Lingji,

I misunderstood your comment. Yes, the extra inputs to the c function are not needed. I have resubmitted the file with these removed for greater speed. Thank you for your suggestion.

Roel H,
Agreed on both counts. The code has been updated and re-posted. Doing some quick testing, the updates you recommended significantly improve speed for very large matrices, thank you.

Zachary,

I'm not quite sure about your question, but let me try to answer anyway:

The "ending point" recursion for the chain should parallel the value recursion which stores the results in CA (without repeating the work already done). In this case you would need a second variable to store where the "last time" ending point was. I am saying "this time" and "last time" because the recursion goes backwards. For example, in the (i > 1 && j > 1) case, the ending point is (i,j), but depending on which of the 3 cases has the min value, the "last time" could be (i-1, j-1), (i-1, j) or (i, j-1).

Hope this helps.

Zachary,

Here are some more suggestions:

1. You can add a second return value that gives the coupling sequence. More on this below.

2. You can move the distance function handle into the call list, as an optional third argument. This makes it more versatile without incurring extra cost.

3. To get the coupling sequence, you can add a cell array at the level of CA (by the way, persistence is not needed), with each cell recording two things: an index into the cell array that gives where the coupling sequence ends "last time," and the new coupling pair that is added "this time."

4. You can change the c() function to do two things: modifying CA, and modifying the cell array, during the if-then-else branches. Return value is not needed.

5. After a call to c(p, q), you can get the distance from CA(p, q), and trace backwards the coupling sequence starting from the cell array at (p, q).

Hope this is clear enough.