Newton-type Method for the Orthonormal INDSCAL Problem in Metric Multidimensional Scaling

The orthonormal individual differences scaling (O-INDSCAL) model is a fundamental tool for metric multidimensional scaling of multiple doubly centered dissimilarity matrices, but efficient algorithms remain scarce when orthonormality and nonnegativity constraints are imposed in medium- to large-scale settings. In this work, we revisit O-INDSCAL from a Riemannian optimization perspective and make two main contributions. First, starting from explicit closed-form expressions of the Riemannian gradient and Hessian, we represent the Hessian as a Kronecker-product-based linear operator on the tangent space and reformulate the Riemannian Newton equation as a symmetric linear system whose dimension coincides with that of the underlying product manifold. This dimension-reduced formulation enables the use of Krylov subspace methods and substantially lowers the per-iteration cost of Newton steps while preserving second-order accuracy. Second, we design a hybrid algorithm that couples a globally convergent Riemannian curvilinear search with Barzilai--Borwein step sizes and a locally quadratically convergent Riemannian Newton refinement, thereby combining robust globalization with fast local convergence. Extensive numerical experiments on synthetic and real O-INDSCAL benchmarks show that the proposed method attains high-accuracy solutions and compares favorably with projected gradient flows, ALS/MPE-type schemes, and generic Manopt-based Riemannian solvers in terms of both efficiency and reliability. MATLAB implementations of all methods used in the numerical experiments are available at Hybrid Newton for OINDSCAL MATLAB Codes.