error doesn't decrease

[codeprof]

To me it seems there is something wrong with the sammon mapping. The error does not decrease and it seems to do not deliver good results if pca is not used.

[tompollard]

Thanks @codeprof, I'll look into this when I get a chance to see if I can work out what's happened. It may be a result of a recently merged pull request.

[codeprof]

At the moment I use the following code (old) with small changes. I do not use the Heassian (commented out) as this seems to result in problems with PCA initialisation. The values in H are "nan", because the result of PCA seems to be complex numbers...

`def sammon(x, n = 2, display = 2, inputdist = 'raw', maxhalves = 20, maxiter = 500, tolfun = 1e-9, init = 'pca'):

import numpy as np 
from scipy.spatial.distance import cdist

"""Perform Sammon mapping on dataset x

y = sammon(x) applies the Sammon nonlinear mapping procedure on
multivariate data x, where each row represents a pattern and each column
represents a feature.  On completion, y contains the corresponding
co-ordinates of each point on the map.  By default, a two-dimensional
map is created.  Note if x contains any duplicated rows, SAMMON will
fail (ungracefully). 

[y,E] = sammon(x) also returns the value of the cost function in E (i.e.
the stress of the mapping).

An N-dimensional output map is generated by y = sammon(x,n) .

A set of optimisation options can be specified using optional
arguments, y = sammon(x,n,[OPTS]):

   maxiter        - maximum number of iterations
   tolfun         - relative tolerance on objective function
   maxhalves      - maximum number of step halvings
   input          - {'raw','distance'} if set to 'distance', X is 
                    interpreted as a matrix of pairwise distances.
   display        - 0 to 2. 0 least verbose, 2 max verbose.
   init           - {'pca', 'random'}

The default options are retrieved by calling sammon(x) with no
parameters.

File        : sammon.py
Date        : 18 April 2014
Authors     : Tom J. Pollard ([email protected])
            : Ported from MATLAB implementation by 
              Gavin C. Cawley and Nicola L. C. Talbot

Description : Simple python implementation of Sammon's non-linear
              mapping algorithm [1].

References  : [1] Sammon, John W. Jr., "A Nonlinear Mapping for Data
              Structure Analysis", IEEE Transactions on Computers,
              vol. C-18, no. 5, pp 401-409, May 1969.

Copyright   : (c) Dr Gavin C. Cawley, November 2007.

This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA

"""

# Create distance matrix unless given by parameters
if inputdist == 'distance':
    D = x
else:
    D = cdist(x,x)

# Remaining initialisation
N = x.shape[0] # hmmm, shape[1]?
scale = 0.5 / D.sum()
D = D + np.eye(N)
Dinv = 1 / D # Returns inf where D = 0.
Dinv[np.isinf(Dinv)] = 0 # Fix by replacing inf with 0 (default Matlab behaviour).
if init == 'pca':
    [UU,DD,_] = np.linalg.svd(x)
    y = UU[:,:n]*DD[:n] 
else:
    y = np.random.normal(0.0,1.0,[N,n])
one = np.ones([N,n])
d = cdist(y,y) + np.eye(N)
dinv = 1. / d # Returns inf where d = 0. 
dinv[np.isinf(dinv)] = 0 # Fix by replacing inf with 0 (default Matlab behaviour).
delta = D-d 
E = ((delta**2)*Dinv).sum() 

# Get on with it
for i in range(maxiter):

    # Compute gradient, Hessian and search direction (note it is actually
    # 1/4 of the gradient and Hessian, but the step size is just the ratio
    # of the gradient and the diagonal of the Hessian so it doesn't
    # matter).
    delta = dinv - Dinv
    deltaone = np.dot(delta,one)
    g = np.dot(delta,y) - (y * deltaone)
    dinv3 = dinv ** 3
    y2 = y ** 2
    H = np.dot(dinv3,y2) - deltaone - np.dot(2,y) * np.dot(dinv3,y) + y2 * np.dot(dinv3,one)       
    s = -g.flatten(order='F') #/ np.abs(H.flatten(order='F'))
    y_old    = y

    # Use step-halving procedure to ensure progress is made
    for j in range(maxhalves):
        s_reshape = s.reshape(2,len(s)/2).T
        y = y_old + s_reshape
        d = cdist(y, y) + np.eye(N)
        dinv = 1 / d # Returns inf where D = 0. 
        dinv[np.isinf(dinv)] = 0 # Fix by replacing inf with 0 (default Matlab behaviour).
        delta = D - d
        E_new = ((delta**2)*Dinv).sum()
        if E_new < E:
            break
        else:
            s = np.dot(0.5,s)

    # Bomb out if too many halving steps are required
    if j == maxhalves:
        print('Warning: maxhalves exceeded. Sammon mapping may not converge...')

    # Evaluate termination criterion
    if np.abs((E - E_new) / E) < tolfun:
        if display:
            print('TolFun exceeded: Optimisation terminated')
        break

    # Report progress
    E = E_new
    if display > 1:
        print('epoch = ' + str(i) + ': E = ' + str(E * scale))

# Fiddle stress to match the original Sammon paper
E = E * scale

return [y,E]

`

[BenWutzke]

Unfortunately this is common with Quasi-Newton. Since the Hessian is not a true Hessian but an approximation, there is a tendency to bounce around a local or global minimum. Sammon NLM only works well if you have a good initial approximation (as from PCA or boiler MDS schemes). Try keeping track of a few iterations of E and stop if you are seeing them again and again.