Finite Metric Embeddings: Good Algorithm?

I have a finite metric space given as a (symmetric) k by k distance matrix. I would like an algorithm to (approximately) isometrically embed this in euclidean space R^(k-1). While it is not always possible to do exactly by solving the system of equations given by the distances I am looking for a solution that embeds with some (very small) controllable error.

I currently use multidimensional scaling (MDS) with the output dimension set to (k-1). It occurs to me that in general MDS may be optimized for the situation where you are trying to r开发者_如何学JAVAeduce the ambient embedding dimension to something less then (k-1) (typically 2 or 3) and that there may be a better algorithm for my restricted case.

Question: What is a good/fast algorithm for realizing a metric space of size k in R^{k-1} using euclidean distance?

Some parameters and pointers:

(1) My k's are relatively small. Say 3 < k < 25

(2) I don't actually care if I embed in R^{k-1}. If it simplifies things/makes things faster any R^N would also be fine as long as it's isometric. I'm happy if there's a faster algorithm or one with less error if I increase to R^k or R^(2k+1).

(3) If you can point to a python implementation I'll be even happier.

(4) Anything better then MDS would work.

Why not try LLE , you can also find the code there.

Ok, as promised here a simple solution:

Notation: Let d_{i,j}=d_{j,i} denote the squared distance between points i and j. Let N be the number of points. Let p_i denote the point i and p_{i,k} the k-th coordinate of the point.

Let's hope I derive the algorithm correctely now. There will be some explanations afterwards so you can check the derivation (I hate it when many indexes appear).

The algorithm uses dynamic programming to arrive at the correct solution

Initialization:
p_{i,k} := 0 for all i and k

Calculation:
for i = 2 to N do
   sum = 0
   for j = 2 to i - 1 do
     accu = d_{i,j} - d_{i,1} - p_{j,j-1}^2
     for k = 1 to j-1 do
       accu = accu + 2 p_{i,k}*p_{j,k} - p_{j,k}^2
     done
     p_{i,j} = accu / ( 2 p_{j,j-1} )
     sum = sum + p_{i,j}^2
   done
   p_{i,i-1} = sqrt( d_{i,0} - sum )
done

If I have not done any grave index mistakes (I usually do) this should do the trick.

The idea behind this:

We set the first point arbitrary at the origin to make our life easier. Not that for a point p_i we never set a coordinate k when k > i-1. I.e. for the second point we only set the first coordinate, for the third point we only set first and second coordinate etc.

Now let's assume we have coordinates for all points p_{k'} where k'<i and we want to calculate the coordinates for p_{i} so that all d_{i,k'} are satisfied (we do not care yet about any constraints for points with k>k'). If we write down the set of equations we have

d_{i,j} = \sum_{k=1}^N (p_{i,k} - p_{j,k} )^2

Because both p_{i,k} and p_{j,k} are equal to zero for k>k' we can reduce this to:

d_{i,j} = \sum_{k=1}^k' (p_{i,k} - p_{j,k} )^2

Also note that by the loop invariant all p_{j,k} will be zero when k>j-1. So we split this equation:

d_{i,j} = \sum_{k=1}^{j-1} (p_{i,k} - p_{j,k} )^2 + \sum_{k=j}^{k'} p_{i,j}^2

For the first equation we just get:

d_{i,1} = \sum_{k=1}^{j-1} (p_{i,k} - p_{j,k} )^2 + \sum_{k=j}^{k'} p_i{i,1}^2

This one will need some special treatment later.

Now if we solve all binomials in the normal equation we get:

d_{i,j} = \sum_{k=1}^{j-1} p_{i,k}^2 - 2 p_{i,k} p_{j,k} + p_{j,k}^2 + \sum_{k=j}^{k'} p_{i,j}^2

subtract the first equation from this and you are left with:

d_{i,j} - d_{i,1} = \sum_{k=1}^{j-1} p_{j,k}^2 - 2 p_{i,k} p_{j,k}

for all j > 1.

If you look at this you'll note that all squares of coordinates of p_i are gone and the only squares we need are already known. This is a set of linear equations that can easily be solved using methods from linear algebra. Actually there is one more special thing about this set of equations: The equations already are in triangular form, so you only need the final step of propagating the solutions. For the final step we are left with one single quadratic equation that we can just solve by taking one square root.

I hope you can follow my reasoning. It's a bit late and my head is a bit spinning from those indexes.

EDIT: Yes, there were indexing mistakes. Fixed. I'll try to implement this in python when I have time in order to test it.