A Tutorial on Structure Optimized Proximity Scaling (STOPS)

Stress Models

The MDS models we list below can be fitted with the smacofx package.

The first popular type of PS supported in stops is based on the loss function type stress (Kruskall 1964), employing a quadratic loss function.

A general formulation of a loss function based on a quadratic loss for \(i,j = 1, \dots, N\) is \[\begin{equation} \label{eq:stress} \sigma_{PS}(X)=\sigma_{stress}(X)=\sum_{i<j} w^*_{ij}\left[d^*_{ij}(X)-\delta^*_{ij}\right]^2=\sum^N_{i=1}\sum^N_{j=1} h(w_{ij})\left[g\left(d_{ij}(X)\right)-f(\delta_{ij})\right]^2 \end{equation}\]

Here, the \(w_{ij}\) are finite weights, with \(w_{ij}=0\) if the entry is missing and \(w_{ij}\neq 0\) otherwise. There are a number of optimization techniques one can use to solve this optimization problem.

The distances used is usually the Euclidian distance as the distance fitted to the points in the configuration, \[\begin{equation} \label{eq:dist} d_{ij}(X) = ||x_{i}-x_{j}||_2=\left( \sum_{m=1}^M |x_{im}-x_{jm}|^2 \right)^{1/2} \ i,j = 1, \dots, N. \end{equation}\] Typically, the norm used is the Euclidean norm, so \(p=2\). In standard MDS \(g(\cdot)=f(\cdot)=I(\cdot)\), the identity function.

This formulation enables one to express a large number of PS methods many of which are implemented in stops. In stops we allow to use specific choices for \(f(\cdot)\), \(g(\cdot)\) and \(h(\cdot)\) from the family of power transformations so one can fit the following stress models:

• Explicitly normalized stress: \(w^*_{ij}=(\sum_{ij}\delta^{*2}_{ij})^{-1}\), \(\delta_{ij}^*=\delta_{ij}\), \(d^*_{ij}(X)=d_{ij}(X)\)
• Stress-1: \(w^*_{ij}=(\sum_{ij} d^{*2}_{ij}(X))^{-1}\), \(\delta_{ij}^*=\delta_{ij}\), \(d_{ij}(X)^*=d_{ij}(X)\)
• Sammon stress (Sammon 1969): \(w^*_{ij}=\delta^{*-1}_{ij}\) , \(\delta_{ij}^*=\delta_{ij}\), \(d_{ij}(X)^*=d_{ij}(X)\)
• Elastic scaling stress (McGee 1966): \(w_{ij}=\delta^{*-2}_{ij}\), \(\delta_{ij}^*=\delta_{ij}\), \(d_{ij}(X)^*=d_{ij}(X)\)
• S-stress (Takane et al. 1977): \(\delta^*_{ij}=\delta_{ij}^2\) and \(d^*_{ij}(X)=d^2_{ij}(X)\), \(w^*_{ij}=1\)
• R-stress (de Leeuw, 2014): \(\delta^*_{ij}=\delta_{ij}\) and \(d^*_{ij}=d^{2r}_{ij}\), \(w^*_{ij}=1\)
• Power MDS (Buja et al. 2008, Rusch, Mair & Hornik 2021): \(\delta^*_{ij}=\delta_{ij}^\lambda\) and \(d^*_{ij}=d^\kappa_{ij}\), \(w^*_{ij}=1\)
• Power elastic scaling (Rusch, Mair & Hornik 2023): \(w^*_{ij}=\delta^{*-2}_{ij}\), \(\delta^*_{ij}=\delta_{ij}^\lambda\) and \(d^*_{ij}=d^\kappa_{ij}\)
• Power Sammon mapping (Rusch, Mair & Hornik 2023): \(w^*_{ij}=\delta^{*-1}_{ij}\), \(\delta^*_{ij}=\delta_{ij}^\lambda\) and \(d^*_{ij}=d^\kappa_{ij}\)
• Restricted powerstress (Rusch, Mair & Hornik, 2021): \(\delta^*_{ij}=\delta_{ij}^\lambda\), \(d^*_{ij}=d^\kappa_{ij}\) but with the restriction \(\kappa=\lambda\). The \(w^*_{ij}=w_{ij}^\rho\) for arbitrary \(w_{ij}\) (e.g., a function of the \(\delta_{ij}\))
• Power stress (encompassing all previous models; Buja et al. 2008, Rusch, Mair & Hornik 2021): \(\delta^*_{ij}=\delta_{ij}^\lambda\), \(d^*_{ij}=d^\kappa_{ij}\) and \(w^*_{ij}=w_{ij}^\rho\) for arbitrary \(w_{ij}\) (e.g., a function of the \(\delta_{ij}\))
• Approximate power stress (Rusch, Mair & Hornik, 2021): This approximates power stress by \[\begin{equation} \label{eq:apstress} \sigma_{PS}(X)=\sigma_{apstress}(X)=\sum_{i<j} \delta_{ij}^\upsilon \left(\delta_{ij}^\tau - d_{ij}(X)\right)^2 \end{equation}\] where the relation to p-stress is so that \(\upsilon=\rho + 2\lambda(1 -1/\kappa)\) and \(\tau=\lambda/\kappa\). The approximation of p-stress by ap-stress works well in cases when for the \(x_i,x_j\) for which \(w_{ij}\) is large, the error \(\delta_{ij}^\lambda - d_{ij}(X)^\kappa\) is small, so that \(d_{ij}(X)^\kappa\) is approximated reasonably well by \(\delta_{ij}^\lambda\) and, equivalently, \(d_{ij}(X)^\kappa\) can be approximated well by \(d_{ij}(X) \delta_{ij}^{(\lambda (\kappa-1)/\kappa)}\).
• Local MDS (LMDS; Chen & Buja, 2009): Let \(N_k\) define the symmetric set of nearby pairs of points \((i,j)\) so that \((i,j)\in N_k\) iff \(i\) is among the \(k-\)nearest neighbours of \(j\) or the other way round. Then we have a stress model with \[\begin{equation} \delta^*_{ij}=\begin{cases} \delta_{ij} &\mbox{ if } (i,j) \in N_k\\ \sqrt{w}\delta_{\infty} &\mbox{ if } (i,j) \in N_k\\ \end{cases} \end{equation}\] and \[\begin{equation} d^*_{ij}(X)=\begin{cases} d_{ij}(X) &\mbox{ if } (i,j) \in N_k\\ \sqrt{w}d_{ij}(X) &\mbox{ if } (i,j) \in N_k\\ \end{cases} \end{equation}\] where \(\delta_{\infty} \rightarrow \infty\) is a large ``imputed’’ dissimilarity that is constant and \(w\) a small weight. For simplification, one takes \(w \sim 1/\delta_{\infty}\) in the standard lMDS objective, so one uses the tuning parameter \(\tau=2w\delta_{\infty}\) for given \(k\).

Strain Models

The second popular type of PS supported in stops is based on the loss function type . Here the \(\Delta^*\) are a transformation of the \(\Delta\), \(\Delta^*= f (\Delta)\) so that \(f(\cdot)=-(h\circ l)(\cdot)\) where \(l\) is any function and \(h(\cdot)\) is a double centering operation, \(h(\Delta)=\Delta-\Delta_{i.}-\Delta_{.j}+\Delta_{..}\) where \(\Delta_{i.}, \Delta_{.j}, \Delta_{..}\) are matrices consisting of the row, column and grand marginal means respectively. These then get approximated by (functions of) the inner product matrices of \(X\) \[\begin{equation} \label{eq:dist2} d_{ij}(X) = \langle x_{i},x_{j} \rangle \end{equation}\]

We can thus express classical scaling as a special case of the general PS loss with \(d_{ij}(X)\) as an inner product, \(g(\cdot) = I(\cdot)\) and \(f(\cdot)=-(h \circ I)(\cdot)\).

If we again allow power transformations for \(g(\cdot)\) and \(f(\cdot)\) one can fit the following strain models with stops

• Powerstrain (Buja et al. 2008, Rusch, Mair & Hornik 2021): \(\delta^*_{ij}=-h(\delta_{ij}^\lambda)\), \(d^*_{ij}=d_{ij}\) and \(w_{ij}=w_{ij}^\rho\) for arbitrary \(w_{ij}\)

Energy Models

The third type of PS supported in stops is based on the energy type losses of pairwise attraction and repulsion between objects (Chen & Buja, 2014, Rusch, Mair, Hornik 2023a). It is related to graph drawing: if vertices are equated with objects, the edge weights are the \(\delta^*_{ij}\), the node positions in the layout is the \(X\) and the distance between nodes in the layout are \(d_{ij}(X)\).

One can express badness-of-fit as comprising a repulsion part \(\propto -\delta^*_{ij}d^*_{ij}(X)\) and an attraction part \(\propto d^{*\mu}_{ij}(X)\) (Chen & Buja, 2014). This is
\[\begin{equation} \label{eq:energy} \sigma_{PS}(X)=\sigma_{EM}(X) = \sum_{i<j} w^*_{ij} \left(a d^{*\mu}_{ij}(X) - b \delta^*_{ij}d^*_{ij}(X) \right) \end{equation}\] with \(a,b\) being constants and \(d^*_{ij}(X)=g(d_{ij}(X))\) and \(\delta^*_{ij}=f(d_{ij}(X))\). For \(\mu=2, a=1\) and \(b=2\) this is a stress model where terms depending solely on \(\delta^*_{ij}\) are disregarded for finding \(X\).

Specific instances of energy models that can be fitted in stops are:

• Box-Cox MDS (Chen & Buja, 2014). This is a three-parameter family \(\theta=(\mu^{(BC)},\lambda^{(BC)},\rho^{(BC)})^\top\): \[\begin{equation} \label{eq:bcstress} \sigma_{PS}(X)=\sigma_{BCstress}(X) = \sum_{i<j} \delta_{ij}^{\rho} \left( BC_{\mu+\lambda}(d_{ij}(X)) - \delta_{ij}^{\lambda} BC_{\mu}(d_{ij}(X)) \right) \end{equation}\] with \(\mu, \rho\in \mathbb{R}\) and \(\lambda \in \mathbb{R}_+\). Here \(BC_\alpha\) is the one-parameter Box-Cox transformation with parameter \(\alpha\), \[\begin{equation} \label{eq:bc} BC_{\alpha}(d)=\begin{cases} \frac{d^{\alpha}-1}{\alpha} \mbox{ if } \alpha \neq 0\\ \log(d) \mbox{ if } \alpha = 0\\ \end{cases} \end{equation}\] For \(w_{ij}=\delta_{ij}, \mu=\lambda=1\) this is very similar to the stress-equivalent energy MDS model. Note that here the distance transformation functions in the attraction and repulsion parts need not be equal. \(\sigma_{\text{BC}}(X)\) shares close similarities with stress with power transformations (power stress): They will often lead to a similar configurations for the same parameters (even equivalent if the exponents for repulsion and attraction parts are the same).

Model

Following Rusch, Mair & Hornik (2023a) the general idea is that from given observations \(\Delta\) we look for a configuration \(X\). The \(X\) has properties with regards to its structural appearance, which we call c-structuredness for configuration-structuredness. There are different types of c-structuredness people might be interested in (say, how clustered the result is, or that dimensions are orthogonal or if there is some submanifold that the data live on). We developed indices for these types of c-structuredness that capture that essence in the configuration (see Rusch, Mair & Hornik 2023b, or Rusch, Mair & Hornik 2020 for a large list of structures and indices).

We have as part of a STOPS model a proximity scaling loss function \(\sigma_{PS}(\cdot)\) and some transformation \(f_{ij}(\delta_{ij}|\theta_f)\) and \(g_{ij}(d_{ij}|\theta_g)\) and possibly \(h_{ij}(w_{ij}|\theta_h)\). These functions are parametrized and all parameters are collected in the vector \(\theta=(\theta_f,\theta_g,\theta_h)\). \(\theta\) is finite, e.g., a transformation applied to all observations like a power transformation. For example for power stress \(\theta=(\kappa,\lambda,\rho)\) and for Box Cox MDS it is \(\theta=(\mu,\lambda,\rho)\). These transformations achieve a sort of push towards more structure, so different values for \(\theta\) will in general lead to different c-structuredness.

Let us assume we are interested in \(L\) different structural qualities of \(X\) and that we have \(L\) corresponding univariate c-structuredness indices \(I_l(X\vert \gamma)\) for the \(l=1,\dots, L\) different structures, capturing the essence of the structural appearance of the configuration with respect to the \(l\)-th structure. For example, we might be interested in both the structural appearance of how clustered the configuration is (structure 1) and how strongly linearly related the column vectors of the configuration are (structure 2). We then measure the c-structuredness of \(X\) for the two structures with an index for clusteredness and one for linear dependence respectively. The \(\gamma\) are optional metaparameters for the indices, which we assume are given and fixed; they control how c-structuredness is measured. We further assume broadly that the transformations \(f(\Delta\vert\theta_f)\) and/or \(g(D(X)\vert \theta_g)\) and/or \(h(W\vert \theta_h)\) produce different c-structuredness in \(X\) for different values of \(\theta\).

In a nutshell our proposal is to select optimal hyperparameters \(\theta^\ast\) for the scaling procedure by assessing the c-structuredness of an optimal configuration \(X^\ast\) found from a PS method for given \(\theta\) usually in combination with its badness-of-fit value. We aim at finding a \(\theta^\ast\) that, when used as transformation parameters in the PS method, will give a configuration that has high (or low) values of the c-structuredness indices. We view this as a multi-objective optimization problem, where we want to maximize/minimize different criteria (either badness-of-fit, or c-structuredness, or both) over \(\theta\). C-structuredness may this way be induced at a possible expense of fit but we control the expense amount.

To formalize this we explicitly write the building blocks of the objective function used for hyperparameter tuning via STOPS as a function of \(\theta\): Let us denote by \(X^\ast(\theta)\) the optimal solution from minimizing a badness-of-fit \(\sigma_{PS}(X \vert \theta)\) for a fixed \(\theta\), so \(X^\ast(\theta):= \arg\min_{X} \sigma_{PS}(X\vert \theta)\). Further we also have the \(L\) different univariate indices with possible metaparameters \(\gamma\), \(I_l(X^\ast(\theta)\vert \gamma)\) to be optimized for.

Specific variants of STOPS can be instantiated by defining objective functions \(\text{Stoploss}(\theta\vert v_0, \dots, v_L, \gamma)\), comprising either badness-of-fit or c-structuredness indices or both in a scalarized combination. Two variants of objective functions are of the following form:

• Additive STOPS (aSTOPS) \[\begin{equation} \label{eq:astops} \text{Stoploss}_{aSTOPS}(\theta\vert v_0, \dots, v_L, \gamma) = v_0 \cdot \sigma_{PS}(X^\ast(\theta)\vert\theta) + \sum^L_{l=1} v_l I_l(X^\ast(\theta)\vert\gamma) \end{equation}\]

• Multiplicative STOPS (mSTOPS) \[\begin{equation} \label{eq:mstops} \text{Stoploss}_{mSTOPS}(\theta\vert v_0, \dots, v_L, \gamma) = \sigma_{PS}(X^\ast(\theta)\vert\theta)^{v_0} \cdot \prod^L_{l=1} I_l(X^\ast(\theta)\vert\gamma)^{v_l} \end{equation}\]

with \(v_0 \in \mathbb{R}_{\geq 0}\) and \(v_1,\dots,v_L \in \mathbb{R}\) being the scalarization weights that determine how the individual parts (MDS loss and c-structuredness indices) are aggregated. Note that in this formulation the aggregation is a sum/product, so the weights must be negative if a higher index stands for more structure and we want more of that structure. Alternatively, the stressweight can be negative and the strucweight positive.

Numerically, the badness-of-fit function value \(\sigma_{PS}(X^\ast(\theta)\vert \theta)\) needs to be normalized to be scale-free and commensurable for comparability of different values of \(\theta\). The objective function for aSTOPS is fully compensatory, whereas for mSTOPS it ensures that a normalized badness-of-fit of \(0\) will always lead to the minimal \(\text{Stoploss}(\theta\vert v_0, \dots, v_L, \gamma)\) for a positive value of \(I_l(\cdot)\). For notational convenience, we will refer to the objective functions for STOPS variants by \(\text{Stoploss}(\theta)\) from now on.

The job is then to find \[\begin{equation} \arg\min_{\theta}\ \text{Stoploss}(\theta) \end{equation}\]

Minimizing stoploss can be quite difficult. In stops we use a nested algorithm combining optimization that internally first solves for \(X\) given \(\theta\), \(\arg\min_X \sigma_{PS}\left(X|\theta\right)\) in an inner loop, and then optimize over \(\theta\) with a metaheuristic in an outer loop. The number of iterations of the outer loop can be controlled with itmax, the number of iterations of the inner loop with itmaxps.

data(kinshipdelta)
dis<-kinshipdelta
X<-mds(dis)$conf
c_stringiness(X)

## [1] 0.1714678

c_clumpiness(X)

## [1] 0.1857798

c_association(X,alpha=0.7,C=10)

## [1] 0.3982779

c_clusteredness(X,minpts=3)

## [1] 0.1807154

structures<-c("cclusteredness","clinearity","cdependence")

strucpars<-list(
                list(epsilon=10,minpts=2,dmax=1.3),  # element 1: list of arguments for c-clusteredness 
                NULL,                                # element 2: argument for c-linearity (empty as it has no parameters)
                list(index=1.2)                      # element 3: list of arguments for c-dependence
                ) 

set.seed(666)
ressm<-stops(kinshipdelta,loss="stress",stressweight=1,structures=c("cclusteredness","clinearity","cdependence"),strucweight=c(-0.33,0.33,-0.33),verbose=0,strucpars=strucpars,stoptype="multiplicative",lower=0,upper=8)
ressm

## 
## Call: stops(dis = kinshipdelta, loss = "stress", structures = c("cclusteredness", 
##     "clinearity", "cdependence"), stressweight = 1, strucweight = c(-0.33, 
##     0.33, -0.33), strucpars = strucpars, lower = 0, upper = 8, 
##     verbose = 0, stoptype = "multiplicative")
## 
## Model: multiplicative STOPS with ratio stress loss function and theta parameter vector (lambda)  =  1.432493 
## 
## Number of objects: 15 
## MDS loss value: 0.06970603 
## C-Structuredness Indices: cclusteredness 0.24108469 clinearity 0.02696876 cdependence 0.20524114 
## Structure optimized loss (stoploss): 0.05705468 
## MDS loss weight: 1 c-structuredness weights: -0.33 0.33 -0.33 
## Number of iterations of ALJ optimization: 28

plot(ressm)

ressa<-stops(kinshipdelta,loss="stress",stressweight=1,structures=c("cclusteredness","clinearity","cdependence"),strucweight=c(-0.33,0.33,-0.33),verbose=0,strucpars=strucpars,stoptype="additive",lower=0,upper=8)
ressa

## 
## Call: stops(dis = kinshipdelta, loss = "stress", structures = c("cclusteredness", 
##     "clinearity", "cdependence"), stressweight = 1, strucweight = c(-0.33, 
##     0.33, -0.33), strucpars = strucpars, lower = 0, upper = 8, 
##     verbose = 0, stoptype = "additive")
## 
## Model: additive STOPS with ratio stress loss function and theta parameter vector (lambda)  =  1.483125 
## 
## Number of objects: 15 
## MDS loss value: 0.07037108 
## C-Structuredness Indices: cclusteredness 0.24482385 clinearity 0.02667558 cdependence 0.20393552 
## Structure optimized loss (stoploss): -0.06891657 
## MDS loss weight: 1 c-structuredness weights: -0.33 0.33 -0.33 
## Number of iterations of ALJ optimization: 25

plot(ressa)

set.seed(666)
resa<-stops(kinshipdelta,structures=c("cassociation","cmanifoldness","clinearity"),loss="powermds",theta=c(1,1),strucpars=list(NULL,NULL,NULL),strucweight=c(-0.5,-0.5,0.5),lower=c(0,0),upper=c(5,5),optimmethod="tgp",itmaxps=1000,itmax=10)

resa

## 
## Call: stops(dis = kinshipdelta, loss = "powermds", theta = c(1, 1), 
##     structures = c("cassociation", "cmanifoldness", "clinearity"), 
##     strucweight = c(-0.5, -0.5, 0.5), strucpars = list(NULL, 
##         NULL, NULL), optimmethod = "tgp", lower = c(0, 0), upper = c(5, 
##         5), itmax = 10, itmaxps = 1000)
## 
## Model: additive STOPS with ratio powermds loss function and theta parameter vector (kappa lambda)  =  0.4893293 4.298529 
## 
## Number of objects: 15 
## MDS loss value: 0.2369349 
## C-Structuredness Indices: cassociation 0.515506235 cmanifoldness 0.953516033 clinearity 0.002089034 
## Structure optimized loss (stoploss): -0.4965317 
## MDS loss weight: 1 c-structuredness weights: -0.5 -0.5 0.5 
## Number of iterations of tgp optimization: 10

plot(resa)

set.seed(666)
resa2<-stops(kinshipdelta,structures=c("cclusteredness"),loss="powermds",strucpars=list(list(epsilon=10,dmax=NULL,minpts=2)),strucweight=-1,stressweight=0,lower=c(0.1,0.1),upper=c(5,5),theta=c(1,1), optimmethod="tgp",itmax=20,verbose=3)

## Starting Optimization 
##  EGO (TGP) Optimization 
## Evaluating function at initial points 
## Starting Bayesian Optimization 
## Found minimum after 20  iterations at 0.2771 1.4104 with stoploss= -0.6696 and default scaling loss= 0.0612 and c-structuredness indices: cclusteredness 0.669581 . Thanks for your patience.

resa2

## 
## Call: stops(dis = kinshipdelta, loss = "powermds", theta = c(1, 1), 
##     structures = c("cclusteredness"), stressweight = 0, strucweight = -1, 
##     strucpars = list(list(epsilon = 10, dmax = NULL, minpts = 2)), 
##     optimmethod = "tgp", lower = c(0.1, 0.1), upper = c(5, 5), 
##     verbose = 3, itmax = 20)
## 
## Model: additive STOPS with ratio powermds loss function and theta parameter vector (kappa lambda)  =  0.277098 1.410424 
## 
## Number of objects: 15 
## MDS loss value: 0.06121718 
## C-Structuredness Indices: cclusteredness 0.669581 
## Structure optimized loss (stoploss): -0.669581 
## MDS loss weight: 0 c-structuredness weights: -1 
## Number of iterations of tgp optimization: 20

ressa<-stops(kinshipdelta,theta=1,structures=c("clinearity"),strucweight=0,loss="stress",verbose=0,upper=8,lower=0,optimmethod="pso",itmax=100)

COPS

A special STOPS model is P-COPS (Rusch, Mair & Hornik 2021). Along the lines of STOPS the overall objective function, which we call , is a weighted combination of the \(\theta-\)parametrized loss function, \(\sigma_{PS}\left(X(\theta)|\theta\right)\), and a c-clusteredness measure, the OPTICS cordillera or \(OC(X(\theta);\epsilon,k,q)\) to be optimized as a function of \(\theta\) or \[\begin{equation} \label{eq:spstress} \text{coploss}(\theta) = v_1 \cdot \sigma_{PS}\left(X(\theta) | \theta \right) - v_2 \cdot \text{OC}\left(X(\theta);\epsilon,k,q\right) \end{equation}\]

with \(v_1,v_2 \in \mathbb{R}\) controlling how much weight should be given to the scaling fit measure and the c-clusteredness.

The c-clusteredness index we use is the OPTICS cordillera which measures how clustered a configuration appears. It is based on the OPTICS algorithm that outputs an ordering together with a distance. The OPTICS cordillera is now simply an agregation of that information. Since we know what constitutes a maximally clustered result, we can derive an upper bound and normalize the index to lie between 0 and 1. If it is maximally clustered, the index gets a value of 1,and it gets 0 if all points are equidistant to their nearest neighbours (a matchstick embedding). See the vignette in the COPS package for more.

A Full Blown Analysis

We recreate parts of an example of Rusch, Mair and Hornik (2023) here. The data are 500 handwritten pen digits and we subject them to different flexible MDS methods and select hyperparameters with STOPS based on c-clusteredness and c-manifoldness. We use Sammon mapping with power transformations for the proximities, Box-Cox MDS, lMDS (all three like in the article) and also approximate power stress (different from the article). The data are available in the package via data(Pendigits500).

This is a tough data set for our implementations because it is quite large, so an individual MDS takes a while. When running STOPS we have to fit a large number of MDS, so this can take quite a long time and therefore we also time so that users can see what to expect.

# data setup
data(Pendigits500)
pendss <- Pendigits500
names(pendss)[17] <- "digit"

## setup data and the c-structuredness hyperparameters
dis <- dist(pendss[,1:16]) #only features are used to get the dissimilarity matrix, not the digit label

# parameters for OC (c-clusteredness)
q <- 2
rang <- c(0,0.6) 
minpts <- 5
eps <- 10
scale <- 3

We start with an initial solution (a Sammon mapping) and see how the c-structuredness of that result is

## initial sammon mapping solution
initsam <- sammon(dis)

## Initial stress        : 0.15233
## stress after   1 iters: 0.15079

## c-structuredness of initial solution
c1 <- cordillera::cordillera(initsam$points,minpts=minpts,q=q,epsilon=eps,rang=rang,scale=scale)
ccom1 <- c_manifoldness(initsam$points)
c1

##        raw     normed 
## 0.70519137 0.08331615

ccom1

## [1] 0.593913

So the c-clusteredness was 0.0833162 and c-manifoldness was 0.593913. This is our reference solution to see if we can improve c-structuredness.

Next we set up the structures and the parameters for the structuredness indices

# #c-structuredness metaparameters 
strucpars <- list(list(minpts=minpts,epsilon=eps,rang=rang,scale=scale), #cordillera 
                  list(NULL)  #manifoldness
                  ) 
structures <- c("cclusteredness","cmanifoldness") #structures

We use weights of \(-400\) for c-clusteredness and \(-2.5\) for c-manifoldness.

## c-structuredness weights 
strucweight <- c(-400,-2.5)

Now we have everything to run STOPS. We use Sammon mapping with powers for the dissimilarities, local MDS and Box-Cox MDS. The search space is from \(1\) to \(6\) and we use optimmethod='tgp' (one can set verbose=4 to see details of the iteration progress). We also time the calculations.

#Sammon with power transformations of proximities
set.seed(666)
timesam <- system.time(restgp <- stops(dis,theta=5.6,loss="sammon",structures=structures,stressweight=1,strucweight=strucweight,strucpars=strucpars,lower=1,upper=6,optimmethod="tgp",model="btgpllm",itmax=10))

#Box Cox MDS
set.seed(667)
timebc <- system.time(restgpbc <- stops(dis,theta=c(4.63,2.77,5.12),loss="bcstress",structures=structures,stressweight=1,strucweight=strucweight,strucpars=strucpars,lower=c(0.1,0.1,0.1),upper=c(6,6,6),optimmethod="tgp",model="btgpllm",itmax=10,itmaxps=1e4))

#lMDS
set.seed(669)
timelmds <- system.time(restgplmds <- stops(dis,theta=c(14.5,3),loss="lmds",weightmat=dis,structures=structures,stressweight=1,strucweight=strucweight,strucpars=strucpars,lower=c(5,0),upper=c(50,3),optimmethod="tgp",model="btgpllm",itmax=10,itmaxps=1e4))

#approx p stress
set.seed(670)
timeaps <- system.time(restgpaps <- stops(dis,theta=c(1,1,1),loss="apstress",structures=structures,stressweight=1,strucweight=strucweight,strucpars=strucpars,lower=c(0.1,0.1,0.1),upper=c(3,3,3),optimmethod="ALJ",itmax=20,itmaxps=1e4,verbose=4))

We now plot the results. We color the digits differently and adjust all the configurations to the same plotting region via Procrustes transformation (as MDS results are invariant to translation, scaling and rotation).

#plot the digit results
colsopt <- cols1 <- colstraj <- factor(pendss[,17])
levels(colsopt) <- hcl.colors(10,"Dark 3")

#Procrustes adjustment to lmds configuration for STOPS optimized
samconf <- restgp$fit$conf
lmdsconf <- restgplmds$fit$conf
bcconf <- restgpbc$fit$conf
apsconf<-restgpaps$fit$conf
samconf <- Procrustes(lmdsconf,samconf)$Yhat
bcconf <- Procrustes(lmdsconf,bcconf)$Yhat
apsconf <- Procrustes(lmdsconf,apsconf)$Yhat

And here are the plots.

par(mfrow=c(2,2))
#Sammon
plot(samconf[,1],samconf[,2],col=as.character(colsopt),type="n",xlab="",ylab="",asp=1,axes=TRUE,main="Sammon")
points(samconf[,1],samconf[,2],col=as.character(colsopt),pch=as.character(pendss[,17]))
box()

#lmds
plot(lmdsconf[,1],lmdsconf[,2],col=as.character(colsopt),type="n",xlab="",ylab="",asp=1,main="lMDS",axes=TRUE)
points(lmdsconf[,1],lmdsconf[,2],col=as.character(colsopt),pch=as.character(pendss[,17]))
box()

#bc stress
plot(bcconf[,1],bcconf[,2],col=as.character(colsopt),type="n",xlab="",ylab="",asp=1,axes=TRUE,main="BC-MDS")
points(bcconf[,1],bcconf[,2],col=as.character(colsopt),pch=as.character(pendss[,17]))
box()

#ap stress
plot(apsconf[,1],apsconf[,2],col=as.character(colsopt),type="n",xlab="",ylab="",asp=1,axes=TRUE,main="Approx. POST-MDS")
points(apsconf[,1],apsconf[,2],col=as.character(colsopt),pch=as.character(pendss[,17]))
box()

We can see that Sammon mapping with power transformation shows a clusteredness that is both high and sensible. The values is 0.1500428 for Sammon. Box-Cox MDS also give s arelatively clustered result, but it seems it finds 3 clusters and that is less in line wht the ground truth. The value is 0.1302618.

restgp

## 
## Call: stops(dis = dis, loss = "sammon", theta = 5.6, structures = structures, 
##     stressweight = 1, strucweight = strucweight, strucpars = strucpars, 
##     optimmethod = "tgp", lower = 1, upper = 6, itmax = 10, model = "btgpllm")
## 
## Model: additive STOPS with ratio sammon loss function and theta parameter vector (lambda)  =  5.41942 
## 
## Number of objects: 500 
## MDS loss value: 0.4897205 
## C-Structuredness Indices: cclusteredness 0.1500428 cmanifoldness 0.9282179 
## Structure optimized loss (stoploss): -61.84794 
## MDS loss weight: 1 c-structuredness weights: -400 -2.5 
## Number of iterations of tgp optimization: 10

restgpbc

## 
## Call: stops(dis = dis, loss = "bcstress", theta = c(4.63, 2.77, 5.12), 
##     structures = structures, stressweight = 1, strucweight = strucweight, 
##     strucpars = strucpars, optimmethod = "tgp", lower = c(0.1, 
##         0.1, 0.1), upper = c(6, 6, 6), itmax = 10, itmaxps = 10000, 
##     model = "btgpllm")
## 
## Model: additive STOPS with ratio bcstress loss function and theta parameter vector (mu lambda rho)  =  4.603888 2.711544 5.411521 
## 
## Number of objects: 500 
## MDS loss value: 0.08998973 
## C-Structuredness Indices: cclusteredness 0.1302618 cmanifoldness 0.9773570 
## Structure optimized loss (stoploss): -54.45813 
## MDS loss weight: 1 c-structuredness weights: -400 -2.5 
## Number of iterations of tgp optimization: 10

Local MDS and approximate power stress give configurations that are less clustered, although the arrangement of points is still sensible even if they don’t align in clusters.

restgpbc

## 
## Call: stops(dis = dis, loss = "bcstress", theta = c(4.63, 2.77, 5.12), 
##     structures = structures, stressweight = 1, strucweight = strucweight, 
##     strucpars = strucpars, optimmethod = "tgp", lower = c(0.1, 
##         0.1, 0.1), upper = c(6, 6, 6), itmax = 10, itmaxps = 10000, 
##     model = "btgpllm")
## 
## Model: additive STOPS with ratio bcstress loss function and theta parameter vector (mu lambda rho)  =  4.603888 2.711544 5.411521 
## 
## Number of objects: 500 
## MDS loss value: 0.08998973 
## C-Structuredness Indices: cclusteredness 0.1302618 cmanifoldness 0.9773570 
## Structure optimized loss (stoploss): -54.45813 
## MDS loss weight: 1 c-structuredness weights: -400 -2.5 
## Number of iterations of tgp optimization: 10

restgpaps

## 
## Call: stops(dis = dis, loss = "apstress", theta = c(1, 1, 1), structures = structures, 
##     stressweight = 1, strucweight = strucweight, strucpars = strucpars, 
##     optimmethod = "ALJ", lower = c(0.1, 0.1, 0.1), upper = c(3, 
##         3, 3), verbose = 4, itmax = 20, itmaxps = 10000)
## 
## Model: additive STOPS with ratio apstress loss function and theta parameter vector (kappa lambda nu)  =  2.137869 1.423048 1.649153 
## 
## Number of objects: 500 
## MDS loss value: 0.04557652 
## C-Structuredness Indices: cclusteredness 0.08438802 cmanifoldness 0.79844724 
## Structure optimized loss (stoploss): -35.70575 
## MDS loss weight: 1 c-structuredness weights: -400 -2.5 
## Number of iterations of ALJ optimization: 19

This brings us to c-manifoldness which is very high for Sammon, Box-Cox MDS and lMDS. Approximate power stress MDS (POST-MDS) shows c-manifoldness that is a bit less bit still relatively high (it is less mainly due ot being more noisy around the imagined manifold). The manifolds are alos easy to see: for Sammon it seems to be reminiscent of a spiral, for approximate power stress it is a circle, for Box-Cox MDS it is a circle as well (or three spherical clusters) and for lMDS a manifold is also apparent although I don’t think it has a name.

The interesting thing here is that the objects are arrnaged along the manifolds in a way that suggests a structure of how handwritten digits are related, e.g., that 0s and 8s are related and that they are closer to 5s than to 1s. There is also a way to get from one digit to others along the manifold, e.g, in approx. POST-MDS one would get from the 6s to over the 9s to the 3s. This may give a hypothesis about a relationship amongst them.

Extras

We also briefly describe some other functions in the package stops.

set.seed(210485)
fwild <- function (x) 10*sin(0.3*x)*sin(1.3*x^2) + 0.00001*x^4 + 0.2*x+80
res2<-ljoptim(50, fwild,lower=-50,upper=50,adaptive=FALSE,accd=1e-16,acc=1e-16)
res2

## $par
## [1] -15.81515
## 
## $value
## [1] 67.46773
## 
## $counts
## function gradient 
##      463       NA 
## 
## $convergence
## [1] 0
## 
## $message
## NULL

plot(fwild, -50, 50, n = 1000, main = "ljoptim() minimising 'wild function'")
points(res2$par,res2$value,col="red",pch=19)

conf_adjust(conf1,conf2)