---
title: "Higher Partials of fStress. Who Needs Them ?"
author: "Jan de Leeuw"
date: "Version 02, August 06, 2017"
output:
  html_document:
    keep_md: yes
    number_sections: yes
    toc: yes
  pdf_document:
    keep_tex: yes
    number_sections: yes
    toc: yes
    toc_depth: 3
fontsize: 12pt
graphics: yes
bibliography: fStress.bib
abstract: We define *fDistances*, which generalize Euclidean distances, squared distances, and log distances. The least squares loss function to fit fDistances to dissimilarity data is *fStress*. We give formulas and R/C code to compute partial derivatives of orders one to four of fStress, relying heavily on the use of Faà di Bruno's chain rule formula for higher derivatives. 
---
<script type="text/x-mathjax-config">
MathJax.Hub.Config({
  TeX: { equationNumbers: { autoNumber: "AMS" } }
});
</script>
---
<style type="text/css">

body{ /* Normal  */
   font-size: 18px;
}
td {  /* Table  */
   font-size: 18px;
}
h1 { /* Header 1 */
 font-size: 28px;
 color: DarkBlue;
}
h2 { /* Header 2 */
 font-size: 22px;
 color: DarkBlue;
}
h3 { /* Header 3 */
 font-size: 18px;
 color: DarkBlue;
}
code.r{ /* Code block */
  font-size: 18px;
}
pre { /* Code block */
  font-size: 18px
}
</style>


```{r function_code, echo = FALSE}
dyn.load("fStress.so")
source("fStress.R")
source("check.R")
```
```{r packages, echo = FALSE}
options (digits = 10) 
suppressPackageStartupMessages (library (captioner, quietly = TRUE))
figure_nums <- captioner (prefix = "Figure")
```
Note: This is a working paper which will be expanded/updated frequently. All suggestions for improvement are welcome. The directory [gifi.stat.ucla.edu/fStress](http://gifi.stat.ucla.edu/fStress) has a pdf version, the bib file, the complete Rmd file with the code chunks,  and the R and C source code.

#Introduction

The multidimensional scaling (MDS) loss function *fStress* (@groenen_deleeuw_mathar_C_95) is defined as
\begin{equation}\label{E:fstress}
\sigma(x):=\mathop{\sum\sum}_{1\leq i<j\leq n}w_{ij}(\delta_{ij}-f(x'A_{ij}x))^2
\end{equation}
for some real-valued $f$. The $w_{ij}$ are positive weights, the $\delta_{ij}$ are *dissimilarities*. This does not necessarily imply that they are actual dissimilarity judgments or measurements, they can be arbitrary transformations of such judgments or measurements as well. Both $w_{ij}$ and $\delta_{ij}$ are fixed and known numbers, the $A_{ij}$ are fixed and known symmetric matrices.  Note there is no constraint that either the dissimilarities $\delta_{ij}$ or the fDistances $f(x'A_{ij}x)$ are non-negative, or that $f$ is increasing. Thus minimizing fStress can also be used, for example, to fit fDistances to similarities.

Our notation is somewhat different from the standard MDS notation, so let's make some translations. We use $x=\mathbf{vec}(X)$, with $X$ the MDS configuration of $n$ points in $p$ dimensions. If we define $\nu(i,s)=(s-1)*n+i$, then
$\{X\}_{is}=x_{\nu(i,s)}$. For most of our formulas it is easier to work with vectors than it is to work with matrices, so we will use
$x$ instead of $X$. It is also helpful for our software development, since matrices and higher-dimensional arrays are stored as
vectors anyway.

The $np\times np$ matrices $A_{ij}$ are defined using unit vectors $e_i$ and $e_j$ of length $n$, all zero except for one element that is equal to one. If you like, the $e_i$ are the columns of the identity matrix. We first define $n\times n$ matrices $\mathcal{A}_{ij}$ by
$$
\mathcal{A}_{ij}:=(e_i-e_j)(e_i-e_j)',
$$
and use $p$ copies of $\mathcal{A}_{ij}$ to make the direct sum
$$
A_{ij}:=\underbrace{\mathcal{A}_{ij}\oplus\cdots\oplus\mathcal{A}_{ij}}_{p\text{ times}}.
$$
Thus the squared Euclidean distance between point $i$ and $j$ in the configuration $X$ is $\mathbf{tr}\ X'\mathcal{A}_{ij}X=x'A_{ij}x$. 

##fDistances

An *fDistance* is a function of the form $g(x)=f(x'Ax)$ for some real-valued $f$. Thus minimzinf fStress is fitting fDistances to dissimilarities, using weighted least squares. fStress properly generalizes the usual *stress* function initially proposed by @kruskal_64a and @kruskal_64b, which uses $g(x)=\sqrt{x'Ax}$. The *sStress* loss function of @takane_young_deleeuw_A_77 uses the identity $g(x)=x'Ax$, and the *lStress* function of @ramsay_77 and @ramsay_82 uses the logarithm $g(x)=\log(x'Ax)$.

For $g$ a general power, i.e. $g(x)=(x'Ax)^r$, we get *rStress*, sometimes also known as *qStress*, studied in a host of (unpublished, my fault) papers by @deleeuw_groenen_pietersz_U_06, @groenen_deleeuw_U_10, @deleeuw_U_14c, @deleeuw_groenen_mair_E_16a, @deleeuw_groenen_mair_E_16b, @deleeuw_groenen_mair_E_16c, @deleeuw_groenen_mair_E_16h. URL's and/or DOI's are in the references section.

In @groenen_deleeuw_mathar_C_95 first and second partials of fStress are given. We add third and fourth order partials. It is unclear if the higher order partials have any practical applications. In @deleeuw_groenen_mair_E_16c there are some applications of the second derivatives of rStress. We briefly mention some possible applications of our higher order partials to majorization (a.k.a. MM) algorithms.

#Partials, Partials Everywhere

##Derivatives of fDistances

We give expressions for the first four derivatives of an arbitrary, but sufficiently many times differentiable, fDistance. In fact, we first address a more general problem, which in its turn is a special case of the first four terms of the multivariate Faà di Bruno formula (see, for instance, @constantine_savits_96, @leipnik_pearce_07). We then apply those results to fDistances.

##Faà di Bruno

Suppose $h:\mathcal{R}^n\rightarrow\mathcal{R}$, $g:\mathcal{R}^n\rightarrow\mathcal{R}$, and $f:\mathcal{R}\rightarrow\mathcal{R}$, such that $h(x)=f(g(x))$. We will give expressions for the partials of $h$ of orders up to four. Our formulas separate the parts that depend on $f$ from the parts that depend only on $g$. 

###First

The first partials can be written in the form
$$
\mathcal{D}_ih(x)=\mathcal{D}f(g(x))h_{11}(x),
$$
with
$$
h_{11}(x):=\mathcal{D}_ig(x).
$$

###Second

Next
$$
\mathcal{D}_{ij}h(x)=\mathcal{D}f(g(x))h_{21}(x)+\mathcal{D}^2f(g(x))h_{22}(x),
$$
\begin{align*}
h_{21}(x)&:=\mathcal{D}_{ij}g(x),\\
h_{22}(x)&:=\mathcal{D}_ig(x)\mathcal{D}_jg(x).
\end{align*}

###Third

Next
$$
\mathcal{D}_{ijk}h(x)=\mathcal{D}f(g(x))h_{31}(x)+\mathcal{D}^2f(g(x))h_{32}(x)+\mathcal{D}^3f(g(x))h_{33}(x),
$$
with
\begin{align*}
h_{31}(x)&:=\mathcal{D}_{ijk}g(x),\\
h_{32}(x)&:=\mathcal{D}_jg(x)\mathcal{D}_{ik}g(x)+\mathcal{D}_ig(x)\mathcal{D}_{jk}g(x)+\mathcal{D}_kg(x)\mathcal{D}_{ij}g(x),\\
h_{33}(x)&:=\mathcal{D}_ig(x)\mathcal{D}_jg(x)\mathcal{D}g_k(x).
\end{align*}

###Fourth

And finally
$$
\mathcal{D}_{ijkl}h(x)=\mathcal{D}f(g(x))h_{41}(x)+\mathcal{D}^2f(g(x))h_{42}(x)+\mathcal{D}^3f(g(x))h_{43}(x)+\mathcal{D}^4f(g(x))h_{44}(x),
$$
with
\begin{align*}
h_{41}(x)&:=\mathcal{D}_{ijkl}g(x),\\
h_{42}(x)&:=\mathcal{D}_{jl}g(x)\mathcal{D}_{ik}g(x)+\mathcal{D}_{il}g(x)\mathcal{D}_{jk}g(x)\\
&+\mathcal{D}_jg(x)\mathcal{D}_{ikl}g(x)+\mathcal{D}_ig(x)\mathcal{D}_{jkl}g(x)+\mathcal{D}_kg(x)\mathcal{D}_{ijl}g(x)+\mathcal{D}_lg(x)\mathcal{D}_{ijk}g(x),\\
h_{43}(x)&:=\mathcal{D}_{il}g(x)\mathcal{D}_jg(x)\mathcal{D}g_k(x)+\mathcal{D}_ig(x)\mathcal{D}_{jl}g(x)\mathcal{D}g_k(x)+\mathcal{D}_ig(x)\mathcal{D}_jg(x)\mathcal{D}_{kl}g(x)\\
&+\mathcal{D}_jg(x)\mathcal{D}_{ik}g(x)\mathcal{D}_lg(x)+\mathcal{D}_ig(x)\mathcal{D}_{jk}g(x)\mathcal{D}_lg(x)+\mathcal{D}_kg(x)\mathcal{D}_{ij}g(x)\mathcal{D}_lg(x),\\
h_{44}(x)&=\mathcal{D}_ig(x)\mathcal{D}_jg(x)\mathcal{D}_kg(x)\mathcal{D}g_l(x).
\end{align*}

###Quadratic Forms

If $g(x)=x'Ax$ for some symmetric $A$, as it is in MDS, then $\mathcal{D}_ig(x)=2\{Ax\}_i$ and $\mathcal{D}_{ij}g(x)=2a_{ij}$. Also $\mathcal{D}_{ijk}g(x)=0$ and $\mathcal{D}_{ijkl}g(x)=0$. Note that we write $\{Ax\}_i$ for element $i$ of the vector $Ax$. Also note that the results in this section apply for any symmetric matrix $A$, even if it is not positive-semidefinite. 

For the first partials
$$
\mathcal{D}_ih(x)=2\mathcal{D}f(x'Ax)\{Ax\}_i,
$$
for the second partials
$$
\mathcal{D}_{ij}h(x)=2\mathcal{D}f(x'Ax)a_{ij}+4\mathcal{D}^2f(x'Ax))\{Ax\}_i\{Ax\}_j,
$$
and for the third partials
$$
\mathcal{D}_{ijk}h(x)=4\mathcal{D}^2f(x'Ax)\{\{Ax\}_ia_{jk}+\{Ax\}_ja_{ik}+\{Ax\}_ka_{ij}\}+8\mathcal{D}^3f(x'Ax)\{Ax\}_i\{Ax\}_j\{Ax\}_k.
$$
For the fourth partials we again write
$$
\mathcal{D}_{ijkl}h(x)=\mathcal{D}f(g(x))h_{41}(x)+\mathcal{D}^2f(g(x))h_{42}(x)+\mathcal{D}^3f(g(x))h_{43}(x)+\mathcal{D}^4f(g(x))h_{44}(x),
$$
and now we have
\begin{align*}
h_{41}(x)&:=0,\\
h_{42}(x)&:=4\{a_{jl}a_{ik}+a_{il}a_{jk}\},\\
h_{43}(x)&:=8\{a_{il}\{Ax\}_j\{Ax\}_k+a_{jl}\{Ax\}_i\{Ax\}_k+a_{kl}\{Ax\}_i\{Ax\}_j+\\
&+a_{ik}\{Ax\}_j\{Ax\}_l+a_{jk}\{Ax\}_i\{Ax\}_l+a_{ij}\{Ax\}_k\{Ax\}_l\},\\
h_{44}(x)&:=16\{Ax\}_i\{Ax\}_j\{Ax\}_k\{Ax\}_l.
\end{align*}


##Derivatives of Stress

If we expand fStress we find, using notation used initially by @deleeuw_C_77,
$$
\sigma(x)=C-\rho(x)+\eta(x),
$$
with
$$
\rho(x):=\mathop{\sum\sum}_{1\leq i<j\leq n}w_{ij}\delta_{ij}f(x'A_{ij}x),
$$
$$
\eta(x):=\frac12\mathop{\sum\sum}_{1\leq i<j\leq n}w_{ij}f^2(x'A_{ij}x),
$$
and with $C$ a constant not depending on $x$.

Clearly both $\rho$ and $\eta^2$ are weighted sums of fDistances, where the fDistances in $\eta$ are the squares of the fDistances
in $\rho$. Thus the partial derivatives are also weighted sums of the partial derivatives of the fDistances. Thus the previous formulas give us all we need.

#Implementation

As in many of our previous recent reports, most of the computing is done in a C dialect conforming with the .C() interface of R,
but also fairly easible integrated into C programs that do not depend on the presence of R.

At the moment we have implemented the following fDistances, which all are arbitrary powers of five base functions.
\begin{align}
h(x)&=\{\log(x'Ax)\}^p,\\
h(x)&=\{x'Ax\}^p,\\
h(x)&=\{\exp(x'Ax)\}^p,\\
h(x)&=\left\{\frac{x'Ax}{1+x'Ax}\right\}^p,\\
h(x)&=\{\log(1+x'Ax)\}^p.
\end{align}
There are helper functions in C which evaluate the base functions and all four of their partials at a given point. There is also
R glue for each of the helper functions, and the R function `checkD()` uses the symbolic differentiation capabilities of R to check the helper functions.

Partial derivatives of fDistances are evaluated using the multivariate Faà di Bruno results in the previous sections for the quadratic case with $g(x)=x'Ax$. The derivatives of the powers of the five base functions are computed from the derivatives of the base functions, using the C function `fStressPower()`, which implements a univariate version of Faà di Bruno. Since our formulas cover arbitrary powers of the five base functions they apply to both the fDistances in $\rho$ and their squares in $\eta$.

The two basic C functions are
`faa_di_bruno()` and `fStressFaaDiBruno()`. The first computes the first four partials of $f(x'Ax)$ for general $A$, the second those
of $f(x'Ax)$ for the $A_{ij}$ used in MDS. In the second case we do not have to store $A$ because there is a simple recipe to compute
its elements whenever they are needed. Our C implementation is far from optimal, because we store full multidimensional arrays of partials without using their super-symmetry. We also fill the arrays with full loops over all indices, without using the sparseness of
the $A_{ij}$.

The R function `checkF()` checks the `faa_di_bruno()` derivatives using numeric differentiation (@gilbert_varadhan_16), and the R function `faaR()` provides the glue for `faa_di_bruno`.

Finally there is the C function `fStressPartials()`, which makes the weighted sum of the partials of the fDistances. It has the R glue function `fStressR()`, and for the gradient and the Hessian it is checked against numerical derivatives by `checkS()`. In our examples the 
first and second partials from `fStressR()` are the same as those from the `numDeriv` package. Since the higher derivatives of the 
base functions and their powers given the same results as those from `deriv()` and `D()`, we can be reasonable sure that the higher derivatives are also correct. Checking this, and possibly improving the efficiency of the code, are topics for further research (and for a future version of this paper).

#Discussion

There is clearly some merit to the idea of fDistances, and to the general form of fStress. It covers many of the previous forms of least squares MDS. The first and second partials are not new. They can be used in gradient, Gauss-Newton, and Newton-Raphson minimization algorithms and to draw pseudo-confidence ellipsoids (@deleeuw_E_17q). The third and fourth order partials allow for a better local approximation of fStress, which at least theoretically can lead to faster algorithms and more precise confidence regions.

In future research we will try to exploit convexity and concavity, and bounds on the derivatives, of fDistances and the implications for majorization algorithms to minimize fStress. 

#Appendix: Code

##R Code

###fStress.R

```{r file_auxilary, code = readLines("fStress.R")}
```

###check.R

```{r file_auxilary, code = readLines("check.R")}
```

##C Code

###fStress.h

```{c file_auxilary, code = readLines("fStress.h"), eval = FALSE}
```


###fStress.c

```{c file_auxilary, code = readLines("fStress.c"), eval = FALSE}
```


References