The classical INdividual Differences SCALing (INDSCAL) model is a standard framework for simultaneous metric multidimensional scaling (MDS) of multiple dissimilarity matrices. Its direct INDSCAL variant works directly with squared dissimilarities on the product manifold \({\mathcal {O}}_{0}(n,r)\times {{\textbf{D}}}(r)^{m}\) and naturally accommodates missing entries via indicator-weighted residuals. While this formulation and its Riemannian structure are well documented in the literature, comparatively little attention has been given to simple, first-order algorithms tailored to this setting with missing data. In this paper, we revisit direct INDSCAL fitting with missing values from a Riemannian optimization perspective and develop a streamlined Riemannian gradient method equipped with a Zhang-Hager-type nonmonotone line search. The proposed scheme avoids both Riemannian Hessian evaluations and vector transport operations, admits global convergence guarantees under standard assumptions, and has low per-iteration cost. Extensive numerical experiments on synthetic datasets with systematically controlled missing rates and on two real data sets show that the proposed method delivers solution quality comparable to that of classical projected gradient flows, several first- and second-order solvers from the Manopt toolbox, and three representative Riemannian conjugate-gradient algorithms, while often achieving lower runtime and stable performance across a range of missing-data patterns. These results indicate that the proposed first-order scheme is a practical and efficient alternative for medium- to large-scale direct INDSCAL problems with incomplete dissimilarity information.
