patchGNU Octave - Patches: patch #8559, statistics package: rewrite of...

 
 

patch #8559: statistics package: rewrite of cmdscale for improved matlab compatibility, mathematical accuracy

Submitter:  JD Walsh <jdwalsh>
Submitted:  Sun 26 Oct 2014 03:27:00 AM UTC
   
 
Category:  None Priority:  5 - Normal
Status:  Done Privacy:  Public
Assigned to:  None Open/Closed:  Closed
* Mandatory Fields

Add a New Comment Rich Markup
   

Sun 02 Nov 2014 06:23:05 PM UTC, comment #5: 

Thanks. I applied the send patch as well. This can be closed.

Arno Onken <asnelt>
Sun 02 Nov 2014 05:07:55 PM UTC, comment #4: 

Quick clarification: since the output `e' is meant to be the raw eigenvalues output by `eig', the fixed `cmdscale' does not change the smallest eigenvalue of `e' to 0. It only changes the set of coordinates given in the output `X'. (Similarly, if `e' includes negative or complex eigenvalues, those are unchanged in the output.)

JD Walsh <jdwalsh>
Sun 02 Nov 2014 04:45:14 PM UTC, comment #3: 

I made additional improvements to `cmdscale' and uploaded a new patch. Since the initial version has already been added to the repository, `cmdscalefix.patch' includes only changes made between the initial version and this one.

The initial version of `cmdscale' I uploaded contains a bug. When given an n by n distance matrix (or a vector corresponding to such a matrix), the coordinates returned by `cmdscale' should be no greater than (n-1)-dimensional. However, because of machine approximation `cmdscale' sometimes returns n-dimensional coordinates. For an example, download the attached file `testdim.mat' and using the initial version of `cmdscale', do:

> load 'testdim.mat'; [X, e] = cmdscale (D);


Viewing the results shows that `e' contains a small positive eigenvalue instead of the 0 eigenvalue it should have, and X is 10x10 instead of 10x9, as it should be.

The initial version of `cmdscale' also needed more input testing and more informative error messages than `print_usage()'.

I added code ensuring that if `cmdscale' finds n positive eigenvalues, it throws out the smallest one, since that eigenvalue corresponds to the approximated zero. This can be tested with `testdim.mat'. I also added more robust input checking and descriptive error messages. Here are some tests to illustrate those:

Wrong data types:

> cmdscale (['a','b'; 'b','c']);
> B{1} = 0; B{2} = 0; cmdscale (B);


Wrong matrix forms:

> cmdscale ([1,2,3; 4,5,6]);
> cmdscale ([0,1; 2,0]);
> cmdscale (zeros(3, 3, 3));


Negative entries:

> cmdscale ([0,-1; -1,0]);


Invalid distance matrix:

> cmdscale ([1,2; 2,1]);


(file #32355, file #32356)

JD Walsh <jdwalsh>
Thu 30 Oct 2014 10:48:06 PM UTC, comment #2: 

Thanks, very nice work. I applied your patch.

Arno Onken <asnelt>
Tue 28 Oct 2014 04:44:55 PM UTC, comment #1: 

Adding maintainer of the statistics package to CC list.

Carnë Draug <carandraug>
Group Member
Sun 26 Oct 2014 03:27:00 AM UTC, original submission:  

I'm a PhD candidate in Mathematics at Georgia Tech, and right now that involves a LOT of Octave programming. I came across a problem with the cmdscale function in the Octave-Forge Statistics package, and I'd like to offer this replacement. As time permits, I hope to help write / improve other parts of Octave.

--JD

TL;DR `cmdscale' needed a complete rewrite. The attached patch is a replacement I wrote.

-----

The current `cmdscale' function does NOT perform classical multidimensional scaling (MDS). Its implementation is nothing like the standard equations for MDS, and the resulting matrix Y is completely different that that returned by MatLab's `cmdscale'. Comparing loss functions such as stress and strain give completely different results when using Octave and MatLab, so the current function does not create equivalent solutions. As a result, code written using MatLab's `cmdscale' gives dramatically different results when run in Octave. MatLab's function also supports non-Euclidean matrices (and difference matrices and similarity matrices). The current `cmdscale' does not.

These issues may be happening because the current `cmdscale' function squares the distance matrix and re-centers it on the first row/column. Centering the distance matrix on an arbitrarily-chosen row/column is not the same as the double centering operation which is key to classical MDS. But even with Euclidean matrices, the current `cmdscale' has a tendency to fail its own test assertion when the given matrix sizes are changed. For instance:

for i = 1:400   % eventually fails, though 50 is generally sufficient
  m=4; n=2; X=rand(m,n); D=pdist(X); assert(pdist(cmdscale(D)), D, m*n*eps);
end

In fairness, I suspect the test failure is due to expecting O(eps) accuracy for matrix operations, since O(sqrt(eps)) accuracy is often the best one can achieve when calling `eig'. Thus, any testing errors may not be relevant to the discrepancies above.

I believe one cause of all these issues is made clear in the current documentation, where it states that classical MDS is ``also known as principal coordinate analysis.'' This is not correct, though the confusion is understandable: classical MDS was called principal coordinate analysis (PCO) at one time (Gower 1966). However, there is a subtle yet important difference between the two techniques: when performing the centering operation in PCO the first principal component (representing a baseline/mean) is usually removed (Heiser & Meulman 1983). This does not take place in the double centering operation used for classical MDS. Among other things, this means the coordinates resulting from classical MDS are invariant up to rotation when permuting the distance matrix, while those of PCO are not. (The coordinates of classical MDS are also nested, and they tend to cluster around their centroid.) These differences certainly explain the discrepancy with the MatLab output, and could possibly be connected to any
issues related to the test assertion.

The replacement I have provided follows the standard MDS implementation (Borg & Groenen 2005), and its results match those of MatLab up to column reflection (i.e., a specific column in the output matrix may equal -1 * the same column in MatLab's implementation). Since MDS is only intended to approximate distance, not direction, my results are indistinguishable from MatLab's up to at least O(sqrt(eps)) (when used appropriately for MDS), and so should be compatible with programs dependent on MatLab's `cmdscale' function.

Classical MDS removes all eigenvalues which are not positive real values. The current `cmdscale' checks for and removes negative eigenvalues, but not complex eigenvalues. Truly symmetric matrices never have complex eigenvalues, but their machine approximations in Octave occasionally generate complex near-zero eigenvalues. I added code to test for and remove these.

I also added error checking (the full form of the input matrix must be symmetric and have all nonnegative entries, with zeros on the main diagonal) and I added the option, present in MatLab's code, to calculate MDS from a similarity matrix (in which case the full form of the input matrix must be symmetric and have all nonnegative entries less than or equal to one, with ones on the main diagonal). Sparse and non-Euclidean matrices work just fine fine.

I greatly expanded and (hopefully) clarified the documentation, and I also changed the reference cited in the documention to something more current and -- for my implementation -- more relevant. The book I cited (Borg & Groenen 2005) describes the formulation of classical MDS and was instrumental in my revision of `cmdscale'.

-----

References:

Borg, I. & Groenen, P.J.F. (2005) Modern Multidimensional Scaling, 2nd edition, Springer.

Gower, J.C. (1966) Some distance properties of latent root and vector methods used in multivariate analysis. Biometrika, 53, 325--338.

Heiser, W.J. & Meulman, J. (1983) Analyzing rectangular tables by joint and constrained multidimensional scaling, J. Econometrics, 22, 139–167.

JD Walsh <jdwalsh>

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:
   
   
Comment:
   

Attached Files
file #32355:  cmdscalefix.patch added by jdwalsh (2KiB - text/x-patch)
file #32356:  testdim.mat added by jdwalsh (2KiB - application/octet-stream)
file #32320:  cmdscale.patch added by jdwalsh (8KiB - text/x-patch)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by asnelt (Posted a comment)
  • -email is unavailable- added by carandraug (Arno Onken - maintainer of the statistics package)
  • -email is unavailable- added by jdwalsh (Submitted the item)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only logged-in users can vote.

     

    Follow 7 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2014-11-02 carandraug StatusNone Done
        Open/ClosedOpen Closed
    2014-11-02 jdwalsh Attached File- Added cmdscalefix.patch, #32355
        Attached File- Added testdim.mat, #32356
    2014-10-28 carandraug Summaryrewrite of cmdscale for improved matlab compatibility, mathematical accuracy statistics package: rewrite of cmdscale for improved matlab compatibility, mathematical accuracy
        Carbon-Copy- Added -email is unavailable-
    2014-10-26 jdwalsh Attached File- Added cmdscale.patch, #32320

    Back to the top

    Powered by Savane 3.13-758e.
    Corresponding source code