Convergence of lsqisotonic (in mdscale)

Hi, I am trying to figure out why lsqisotonic does not converge in mdscale (in both stress and sstress criterions) after many hours on a powerful machine. I changed the code a bit to print the number of non-monotonic blocks, and it seems to get stuck at some point.

Is convergence guaranteed? Is there any approaches to speed-up / run this method in parallel?

I would appreciate if someone could direct me to an explanation of this method in laymen's terms. I am new to MDS, so please be gentle.

Here is the part that gets stuck:

diffs = diff(yhat);
while any(diffs<0) % If all blocks are monotonic, then we're done.
    diffs = diff(yhat);
    % Otherwise, merge blocks of non-increasing fitted values, and set the
    % fitted value within each block equal to a constant, the weighted mean
    % of values in that block.
    idx = cumsum([1; (diffs>0)]);
    sumyhat = accumarray(idx,w.*yhat);
    w = accumarray(idx,w);
    yhat = sumyhat ./ w;
    block = idx(block);
    display(horzcat(num2str(sum(diffs<0)),' non monotonic blocks left. magnitude=',num2str(norm(yhat))))
end

Command window:

...

10 non monotonic blocks left. magnitude=118.9195
10 non monotonic blocks left. magnitude=118.9161
10 non monotonic blocks left. magnitude=118.9127
10 non monotonic blocks left. magnitude=118.9093
10 non monotonic blocks left. magnitude=118.9059
10 non monotonic blocks left. magnitude=118.9026
10 non monotonic blocks left. magnitude=118.8992
10 non monotonic blocks left. magnitude=118.8958
10 non monotonic blocks left. magnitude=118.8924
10 non monotonic blocks left. magnitude=118.8891
10 non monotonic blocks left. magnitude=118.8857
10 non monotonic blocks left. magnitude=118.8823
10 non monotonic blocks left. magnitude=118.8789
10 non monotonic blocks left. magnitude=118.8756
10 non monotonic blocks left. magnitude=118.8722
10 non monotonic blocks left. magnitude=118.8688
10 non monotonic blocks left. magnitude=118.8655
10 non monotonic blocks left. magnitude=118.8621
10 non monotonic blocks left. magnitude=118.8588
10 non monotonic blocks left. magnitude=118.8554
10 non monotonic blocks left. magnitude=118.852

Answer 1

[Peter Perkins]

Shay, it kind of impossible to say what's going on here without knowing what's in the data. At each loop iteration, the lines

idx = ...
...
yhat = ...

ought to be merging the runs of decreasing values, and decreasing the lengths of yhat and block. That seems to not be happening, but you'll have to look at some of those vectors to really figure out why the call to accumarray doesn't seem to be "accumulating". Look at [yhat(:) diffs(:) idx(:)] perhaps. It may simply be that you have some kind of garbage in your data that you're not mentioning.

[Shay]

What does that mean in the context of the conjugate gradient method? What would be garbage in the data? What should I be looking for in the matrix you suggested?

My data is a (very) noisy frequency matrix, n=12071. I know the original data points were sampled in 3d and I am trying to embed them back there.

Thanks for the help

[Peter Perkins]

Shay, now that I look at your output, it's not even clear from what you've said that the loop itself is running for a large number of iterations, or whether it is just being called a large number of times. There's just not enough information here to get very far.

If it is the former, which is how I interpreted your "getting stuck", then what you would hope to see is that at each iteration of that loop, the number of elements where diffs<0 decreases until there are none. If you look at [yhat(:) diffs(:) idx(:)] as I suggested, just before the first call to accumarray, you should see that runs of decreasing yhat values are all assigned the same index, and that the next time around that set of yhat values are all the same -- setting yhat values within one block is what the last four lines of that loop do. Actually, I guess you'd have to look at [0; diffs(:)] to make it the right length.

This algorithm does nothing more than find runs of decreasing values, and replace them with their mean.

But you should make sure that it really is getting stuck in that loop, and that the problem is not that outer optimization is just taking a lot of steps. This is after all a _very_ high dimensional optimization.

[Shay]

Hi Peter,

My loop is running for a very large number of iterations (i validated this by setting break points in the code). I let this run for the weekend, it seems to have progressed to 2-non-monotonic blocks, so i am inclined to let it run for a few more days and see what happens, since as you say, this is a hard problem.

From what i understand by your description of the algorithm, convergence is NOT guaranteed, am i right?

I will keep you posted. Thanks again.

Answer 2

[Yun Wang]

I replaced the part of the code that gets stuck with the following code:

n = length(yhat);
b = 0; bstart = zeros(1,n); bend = zeros(1,n);
for i = 1:n
    b = b + 1;
    yhat(b) = yhat(i);
    w(b) = w(i);
    bstart(b) = i; bend(b) = i;
    while b > 1 && yhat(b) < yhat(b-1)
        yhat(b-1) = (yhat(b-1) * w(b-1) + yhat(b) * w(b)) / (w(b-1) + w(b));
        w(b-1) = w(b-1) + w(b);
        bend(b-1) = bend(b);
        b = b - 1;
    end
end
idx = zeros(1,n);
for i = 1:b
    idx(bstart(i) : bend(i)) = i;
end
block = idx(block);

This code segments merges blocks in a O(n), non-iterative fashion. The original implementation is O(n) in each iteration and therefore O(n^2) in total, which is why it can get stuck.

You can download the updated lsqisotonic function here: https://www.mathworks.com/matlabcentral/fileexchange/64933-lsqisotonic-x-y-w-