I.6.5.4: Squared Distance Scaling

Suppose we want to minimize $$\sigma(X)=\sum_{i=1}^m\sum_{j=1}^m(\delta_{ij}-d_{ij}^2(X))^2,$$ with $d_{ij}^2(X)=(x_i-x_j)'(x_i-x_j)$ squared Euclidean distance. An augmentation algorithm for this problem, modestly called ELEGANT, was designed by De Leeuw [1975]. That paper was never published and the manuscript is probably lost, but the algorithm was described, discussed, and applied by both Takane [1977] and Browne [1987]

We augment $\sigma$ to $$\sigma(X,\eta)= \sum_{i=1}^m\sum_{j=1}^m\sum_{k=1}^m\sum_{\ell=1}^m (\eta_{ijk\ell}-(x_i-x_j)'(x_k-x_\ell))^2,$$ where we require $\eta_{ijij}=\delta_{ij}$ while the other $\eta_{ijk\ell}$ are free. Define $C\mathop{=}\limits^{\Delta}XX'$ and assume that $X$ is column-centered, i.e. $C$ is doubly-centered. The augmentation works because $$\sum_{i=1}^m\sum_{j=1}^m\sum_{k=1}^m\sum_{\ell=1}^m \left((x_i-x_j)'(x_k-x_\ell)\right)^2=4n^2\mathbf{tr}\ C^2.$$ Also $$\sum_{i=1}^m\sum_{j=1}^m\sum_{k=1}^m\sum_{\ell=1}^m\eta_{ijk\ell}(x_i-x_j)'(x_k-x_\ell)=\mathbf{tr}\ UC,$$ where $$u_{ij}\mathop{=}\limits^{\Delta}\sum_{k=1}^m\sum_{\ell=1}^m (\eta_{ikj\ell}-\eta_{i\ell kj}-\eta_{kij\ell}+\eta_{\ell i kj}).$$ Thus we can minimize the augmented loss function over $X$ for fixed $\eta_{ijk\ell}$ by minimizing $\mathbf{tr}\ (\frac{1}{4n^2}U-XX')^2$. This means computing the $p$ largest eigenvalues and corresponding eigenvectors of $U$.

Minimizing the augmented loss over the $\eta_{ijk\ell}$ for fixed $X$ is $$\eta_{ijk\ell}^{(k)}= \begin{cases} \delta_{ij}&\text{ if } (i,j)=(k,\ell),\\ (x_i^{(k)}-x_j^{(k)})'(x_k^{(k)}-x_\ell^{(k)})&\text{ otherwise}. \end{cases}$$ This is enough information to get the ALS algorithm going. The "elegance" so far is reducing a problem involving multivariate quartics to a sequence of eigenvalue problems. It is distinctly unelegant, however, that the computations need four-dimensional arrays. But it turns out these can easily be gotten rid of. We use

$$\sum_{i=1}^m\sum_{j=1}^m\sum_{k=1}^m\sum_{\ell=1}^m\eta_{ijk\ell}^{(k)}(x_i-x_j)'(x_k-x_\ell)= 2\mathbf{tr}\ C\left(2n^2C^{(k)}+ B(X^{(k)})\right),$$ where $$B(X^{(k)})\mathop{=}\limits^{\Delta}\sum_{i=1}^m\sum_{j=1}^m(\delta_{ij}-d_{ij}^2(X^{(k)}))A_{ij}$$ and $A_{ij}\mathop{=}\limits^{\Delta}(e_i-e_j)(e_i-e_j)'$. Thus we find $X^{(k+1)}$ by computing eigenvalues and eigenvectors of $C^{(k)}+\frac{1}{2n^2} B(X^{(k)})$ and no intermediate computation or storage of the $\eta_{ijk\ell}$ is required.