# Higher Partials of fStress. Who Needs Them ? Jan de Leeuw Version 02, August 06, 2017 --- 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 \label{E:fstress} \sigma(x):=\mathop{\sum\sum}_{1\leq i= 1") if (order == 1) D(expr, name) else DD(D(expr, name), name, order - 1) } dd <- DD(expr, "x", order) x <- value return (eval (dd)) } checkF <- function (x, fNumber, pPower, a) { f <- function (x, fNumber, pPower, a) { g <- sum (a * outer (x, x)) return ((fList[[fNumber]](g)) ^ pPower) } h0 <- f(x, fNumber, pPower, a) h1 <- grad (f, x, fNumber = fNumber, pPower = pPower, a = a) h2 <- hessian (f, x, fNumber = fNumber, pPower = pPower, a = a) return (list ( x = x, h0 = h0, h1 = h1, h2 = h2 )) } checkS <- function (x, w, delta, p, fNumber, pPower) { f <- function (x, w, delta, p, fNumber, pPower) { r <- length (x) n <- round (r / p) d <- as.vector(dist(matrix(x, n, p))) ^ 2 e <- fList[[fNumber]](d) ^ pPower sum (w * (delta - e) ^ 2) / 2 } h0 <- f( x, w = w, delta = delta, p = p, fNumber = fNumber, pPower = pPower ) h1 <- grad ( f, x, w = w, delta = delta, p = p, fNumber = fNumber, pPower = pPower ) h2 <- hessian ( f, x, w = w, delta = delta, p = p, fNumber = fNumber, pPower = pPower ) return (list ( x = x, h0 = h0, h1 = h1, h2 = h2 )) }  ##C Code ###fStress.h c #ifndef FSTRESS_H #define FSTRESS_H #include #include #include #include static inline int VINDEX(const int); static inline int MINDEX(const int, const int, const int); static inline int SINDEX(const int, const int, const int); static inline int TINDEX(const int, const int, const int); static inline int AINDEX(const int, const int, const int, const int, const int); static inline int CINDEX(const int, const int, const int, const int, const int, const int, const int); static inline double SQUARE(const double); static inline double THIRD(const double); static inline double FOURTH(const double); static inline double FIFTH(const double); static inline double MAX(const double, const double); static inline double MIN(const double, const double); static inline int IMIN(const int, const int); static inline int IMAX(const int, const int); static inline int IMOD(const int, const int); static inline int ASEEK(const int, const int, const int, const int, const int, const int); static inline int VINDEX(const int i) { return i - 1; } static inline int MINDEX(const int i, const int j, const int n) { return (i - 1) + (j - 1) * n; } static inline int AINDEX(const int i, const int j, const int k, const int n, const int m) { return (i - 1) + (j - 1) * n + (k - 1) * n * m; } static inline int CINDEX(const int i, const int j, const int k, const int l, const int n, const int m, const int r) { return (i - 1) + (j - 1) * n + (k - 1) * n * m + (l - 1) * n * m * r; } static inline int SINDEX(const int i, const int j, const int n) { return ((j - 1) * n) - (j * (j - 1) / 2) + (i - j) - 1; } static inline int TINDEX(const int i, const int j, const int n) { return ((j - 1) * n) - ((j - 1) * (j - 2) / 2) + (i - (j - 1)) - 1; } static inline double SQUARE(const double x) { return x * x; } static inline double THIRD(const double x) { return x * x * x; } static inline double FOURTH(const double x) { return x * x * x * x; } static inline double FIFTH(const double x) { return x * x * x * x * x; } static inline double MAX(const double x, const double y) { return (x > y) ? x : y; } static inline double MIN(const double x, const double y) { return (x < y) ? x : y; } static inline int IMAX(const int x, const int y) { return (x > y) ? x : y; } static inline int IMIN(const int x, const int y) { return (x < y) ? x : y; } static inline int IMOD(const int x, const int y) { return (((x % y) == 0) ? y : (x % y)); } static inline int ASEEK(const int n, const int p, const int u, const int v, const int i, const int j) { int value = 0; if (abs (i - j) < n) { for (int s = 1; s <= p; s++) { int k = (s - 1) * n; if ((i == (u + k)) && (j == (u + k))) { value = 1; } if ((i == (u + k)) && (j == (v + k))) { value = -1; } if ((i == (v + k)) && (j == (u + k))) { value = -1; } if ((i == (v + k)) && (j == (v + k))) { value = 1; } } } return (value); } #endif /* FSTRESS_H */  ###fStress.c c #include "fStress.h" void fStressLog(const double *x, double *f0, double *f1, double *f2, double *f3, double *f4) { double xx = *x; *f0 = log(xx); *f1 = 1 / xx; *f2 = -1 / SQUARE(xx); *f3 = 2 / THIRD(xx); *f4 = -6 / FOURTH(xx); return; } void fStressIdentity(const double *x, double *f0, double *f1, double *f2, double *f3, double *f4) { double xx = *x; *f0 = xx; *f1 = 1.0; *f2 = 0.0; *f3 = 0.0; *f4 = 0.0; return; } void fStressExponent(const double *x, double *f0, double *f1, double *f2, double *f3, double *f4) { double xx = *x; *f0 = exp(xx); *f1 = exp(xx); *f2 = exp(xx); *f3 = exp(xx); *f4 = exp(xx); return; } void fStressBounded(const double *x, double *f0, double *f1, double *f2, double *f3, double *f4) { double xx = *x, xz = 1.0 + xx; *f0 = xx / xz; *f1 = 1.0 / SQUARE(xz); *f2 = -2.0 / THIRD(xz); *f3 = 6.0 / FOURTH(xz); *f4 = -24.0 / FIFTH(xz); return; } void fStressLogPlusOne(const double *x, double *f0, double *f1, double *f2, double *f3, double *f4) { double xx = *x, xp = 1.0 + xx, xz = 1 / xp; *f0 = log(xp); *f1 = xz; *f2 = -SQUARE(xz); *f3 = 2.0 * THIRD(xz); *f4 = -6.0 * FOURTH(xz); } void (*fStressTable[5])(const double *, double *, double *, double *, double *, double *) = {fStressLog, fStressIdentity, fStressExponent, fStressBounded, fStressLogPlusOne}; void fStressPower(const double *x, const int *fNumber, const double *ppar, double *g0, double *g1, double *g2, double *g3, double *g4) { double f0, f1, f2, f3, f4, pp = *ppar; (void)fStressTable[*fNumber](x, &f0, &f1, &f2, &f3, &f4); *g0 = pow(f0, pp); *g1 = pp * f1 * pow(f0, pp - 1.0); *g2 = pp * (pp - 1.0) * pow(f0, pp - 2.0) * SQUARE(f1); *g2 += pp * pow(f0, pp - 1.0) * f2; *g3 = pp * (pp - 1.0) * (pp - 2.0) * pow(f0, pp - 3.0) * THIRD(f1); *g3 += 3.0 * pp * (pp - 1.0) * pow(f0, pp - 2.0) * f1 * f2; *g3 += pp * pow(f0, pp - 1.0) * f3; *g4 = pp * (pp - 1.0) * (pp - 2.0) * (pp - 3.0) * pow(f0, pp - 4.0) * FOURTH(f1); *g4 += 6.0 * pp * (pp - 1.0) * (pp - 2.0) * pow(f0, pp - 3.0) * SQUARE(f1) * f2; *g4 += 4.0 * pp * (pp - 1.0) * pow(f0, pp - 2.0) * (f1 * f3); *g4 += 3.0 * pp * (pp - 1.0) * pow(f0, pp - 2.0) * SQUARE(f2); *g4 += pp * pow(f0, pp - 1.0) * f4; } void fStressFaaDiBruno(const double *x, const int *n, const int *p, const int *u, const int *v, const int *fNumber, const double *par, double *ax, double *gx, double *h0, double *h1, double *h2, double *h3, double *h4) { double f0, f1, f2, f3, f4; int nn = *n, uu = *u, vv = *v, pp = *p, np = nn * pp; *gx = 0.0; for (int i = 1; i <= np; i++) { ax[VINDEX(i)] = 0.0; } for (int s = 1; s <= pp; s++) { double xuv = x[MINDEX(uu, s, nn)] - x[MINDEX(vv, s, nn)]; ax[MINDEX(uu, s, nn)] = xuv; ax[MINDEX(vv, s, nn)] = -xuv; } for (int i = 1; i <= np; i++) { *gx += ax[VINDEX(i)] * x[VINDEX(i)]; } (void)fStressPower(gx, fNumber, par, &f0, &f1, &f2, &f3, &f4); *h0 = f0; for (int i = 1; i <= np; i++) { h1[VINDEX(i)] = 2.0 * f1 * ax[VINDEX(i)]; } for (int i = 1; i <= np; i++) { for (int j = 1; j <= np; j++) { h2[MINDEX(i, j, np)] = 2.0 * f1 * ASEEK(nn, pp, uu, vv, i, j) + 4.0 * f2 * ax[VINDEX(i)] * ax[VINDEX(j)]; } } for (int i = 1; i <= np; i++) { for (int j = 1; j <= np; j++) { for (int k = 1; k <= np; k++) { h3[AINDEX(i, j, k, np, np)] = 4.0 * f2 * ((ax[VINDEX(i)] * ASEEK(nn, pp, uu, vv, j, k)) + (ax[VINDEX(j)] * ASEEK(nn, pp, uu, vv, i, k)) + (ax[VINDEX(k)] * ASEEK(nn, pp, uu, vv, i, j))); h3[AINDEX(i, j, k, nn, nn)] += 8.0 * f3 * ax[VINDEX(i)] * ax[VINDEX(j)] * ax[VINDEX(k)]; } } } for (int i = 1; i <= np; i++) { for (int j = 1; j <= np; j++) { for (int k = 1; k <= np; k++) { for (int l = 1; l <= np; l++) { h4[CINDEX(i, j, k, l, np, np, np)] = 4.0 * f2 * (ASEEK(nn, pp, uu, vv, j, l) * ASEEK(nn, pp, uu, vv, i, k) + ASEEK(nn, pp, uu, vv, i, l) * ASEEK(nn, pp, uu, vv, j, k)); h4[CINDEX(i, j, k, l, nn, nn, nn)] += 16.0 * f4 * ax[VINDEX(i)] * ax[VINDEX(j)] * ax[VINDEX(k)] * ax[VINDEX(l)]; h4[CINDEX(i, j, k, l, nn, nn, nn)] += 8.0 * f3 * ASEEK(nn, pp, uu, vv, i, l) * ax[VINDEX(j)] * ax[VINDEX(k)]; h4[CINDEX(i, j, k, l, nn, nn, nn)] += 8.0 * f3 * ASEEK(nn, pp, uu, vv, j, l) * ax[VINDEX(i)] * ax[VINDEX(k)]; h4[CINDEX(i, j, k, l, nn, nn, nn)] += 8.0 * f3 * ASEEK(nn, pp, uu, vv, k, l) * ax[VINDEX(j)] * ax[VINDEX(j)]; h4[CINDEX(i, j, k, l, nn, nn, nn)] += 8.0 * f3 * ASEEK(nn, pp, uu, vv, i, k) * ax[VINDEX(j)] * ax[VINDEX(l)]; h4[CINDEX(i, j, k, l, nn, nn, nn)] += 8.0 * f3 * ASEEK(nn, pp, uu, vv, j, k) * ax[VINDEX(i)] * ax[VINDEX(l)]; h4[CINDEX(i, j, k, l, nn, nn, nn)] += 8.0 * f3 * ASEEK(nn, pp, uu, vv, i, j) * ax[VINDEX(k)] * ax[VINDEX(l)]; } } } } return; } void faa_di_bruno(const double *x, const int *n, const int *fNumber, const double *par, const double *a, double *ax, double *gx, double *h0, double *h1, double *h2, double *h3, double *h4) { double f0, f1, f2, f3, f4; int nn = *n; *gx = 0.0; for (int i = 1; i <= nn; i++) { ax[VINDEX(i)] = 0.0; for (int j = 1; j <= nn; j++) { ax[VINDEX(i)] += a[MINDEX(i, j, nn)] * x[VINDEX(j)]; } *gx += ax[VINDEX(i)] * x[VINDEX(i)]; } (void)fStressPower(gx, fNumber, par, &f0, &f1, &f2, &f3, &f4); *h0 = f0; for (int i = 1; i <= nn; i++) { h1[VINDEX(i)] = 2.0 * f1 * ax[VINDEX(i)]; } for (int i = 1; i <= nn; i++) { for (int j = 1; j <= nn; j++) { h2[MINDEX(i, j, nn)] = 2.0 * f1 * a[MINDEX(i, j, nn)] + 4.0 * f2 * ax[VINDEX(i)] * ax[VINDEX(j)]; } } for (int i = 1; i <= nn; i++) { for (int j = 1; j <= nn; j++) { for (int k = 1; k <= nn; k++) { h3[AINDEX(i, j, k, nn, nn)] = 4.0 * f2 * ((ax[VINDEX(i)] * a[MINDEX(j, k, nn)]) + (ax[VINDEX(j)] * a[MINDEX(i, k, nn)]) + (ax[VINDEX(k)] * a[MINDEX(i, j, nn)])); h3[AINDEX(i, j, k, nn, nn)] += 8.0 * f3 * ax[VINDEX(i)] * ax[VINDEX(j)] * ax[VINDEX(k)]; } } } for (int i = 1; i <= nn; i++) { for (int j = 1; j <= nn; j++) { for (int k = 1; k <= nn; k++) { for (int l = 1; l <= nn; l++) { h4[CINDEX(i, j, k, l, nn, nn, nn)] = 4.0 * f2 * (a[MINDEX(j, l, nn)] * a[MINDEX(i, k, nn)] + a[MINDEX(i, l, nn)] * a[MINDEX(j, k, nn)]); h4[CINDEX(i, j, k, l, nn, nn, nn)] += 16.0 * f4 * ax[VINDEX(i)] * ax[VINDEX(j)] * ax[VINDEX(k)] * ax[VINDEX(l)]; h4[CINDEX(i, j, k, l, nn, nn, nn)] += 8.0 * f3 * a[MINDEX(i, l, nn)] * ax[VINDEX(j)] * ax[VINDEX(k)]; h4[CINDEX(i, j, k, l, nn, nn, nn)] += 8.0 * f3 * a[MINDEX(j, l, nn)] * ax[VINDEX(i)] * ax[VINDEX(k)]; h4[CINDEX(i, j, k, l, nn, nn, nn)] += 8.0 * f3 * a[MINDEX(k, l, nn)] * ax[VINDEX(j)] * ax[VINDEX(j)]; h4[CINDEX(i, j, k, l, nn, nn, nn)] += 8.0 * f3 * a[MINDEX(i, k, nn)] * ax[VINDEX(j)] * ax[VINDEX(l)]; h4[CINDEX(i, j, k, l, nn, nn, nn)] += 8.0 * f3 * a[MINDEX(j, k, nn)] * ax[VINDEX(i)] * ax[VINDEX(l)]; h4[CINDEX(i, j, k, l, nn, nn, nn)] += 8.0 * f3 * a[MINDEX(i, j, nn)] * ax[VINDEX(k)] * ax[VINDEX(l)]; } } } } return; } void fStressPartials(const double *x, const double *w, const double *delta, const int *n, const int *p, const int *fNumber, const double *par, double *stress, double *d, double *fd, double *f1, double *f2, double *f3, double *f4) { double par2 = 2 * *par, nn = *n, pp = *p, np = nn * pp, gx, hx, h0, g0; double *ax = (double *)calloc((size_t)np, sizeof(double)); double *h1 = (double *)calloc((size_t)np, sizeof(double)); double *h2 = (double *)calloc((size_t)SQUARE(np), sizeof(double)); double *h3 = (double *)calloc((size_t)THIRD(np), sizeof(double)); double *h4 = (double *)calloc((size_t)FOURTH(np), sizeof(double)); double *g1 = (double *)calloc((size_t)np, sizeof(double)); double *g2 = (double *)calloc((size_t)SQUARE(np), sizeof(double)); double *g3 = (double *)calloc((size_t)THIRD(np), sizeof(double)); double *g4 = (double *)calloc((size_t)FOURTH(np), sizeof(double)); *stress = 0.0; for (int j = 1; j <= nn - 1; j++) { for (int i = j + 1; i <= nn; i++) { int k = SINDEX(i, j, nn); (void)fStressFaaDiBruno(x, n, p, &i, &j, fNumber, par, ax, &hx, &h0, h1, h2, h3, h4); d[k] = hx; fd[k] = h0; *stress += w[k] * SQUARE(delta[k] - fd[k]); (void)fStressFaaDiBruno(x, n, p, &i, &j, fNumber, &par2, ax, &gx, &g0, g1, g2, g3, g4); for (int r = 1; r <= np; r++) { int ind = VINDEX(r); f1[ind] += w[k] * (0.5 * g1[ind] - delta[k] * h1[ind]); for (int s = 1; s <= np; s++) { int ind = MINDEX(r, s, np); f2[ind] += w[k] * (0.5 * g2[ind] - delta[k] * h2[ind]); for (int t = 1; t <= np; t++) { int ind = AINDEX(r, s, t, np, np); f3[ind] += w[k] * (0.5 * g3[ind] - delta[k] * h3[ind]); for (int u = 1; u <= np; u++) { int ind = CINDEX(r, s, t, u, np, np, np); f4[ind] += w[k] * (0.5 * g4[ind] - delta[k] * h4[ind]); } } } } } } *stress /= 2.0; free(ax); free(h1); free(h2); free(h3); free(h4); free(g1); free(g2); free(g3); free(g4); } ` References