Sample-Efficient Geometry Reconstruction from Euclidean Distances using Non-Convex Optimization

Authors: ipsita  ghosh
abiy  tasissa
christian  kümmerle
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-15 21:26:29

Abstract: The problem of finding suitable point embedding or geometric configurations given only Euclidean distance information of point pairs arises both as a core task and as a sub-problem in a variety of machine learning applications. In this paper, we aim to solve this problem given a minimal number of distance samples. 
To this end, we leverage continuous and non-convex rank minimization formulations of the problem and establish a local convergence  
guarantee for a variant of iteratively reweighted least squares (IRLS), which applies if a minimal random set of observed distances is provided.
 As a technical tool, we establish a restricted isometry property (RIP) restricted to a tangent space of the manifold of symmetric rank-$r$ matrices given random Euclidean distance  measurements, which might be of independent interest for the analysis of other non-convex approaches.  Furthermore, we assess data efficiency, scalability and generalizability of different reconstruction algorithms through numerical experiments with simulated data as well as real-world data, demonstrating the proposed algorithm's ability to identify the underlying geometry from fewer distance samples compared to the state-of-the-art.
 
 The Matlab code can be found at \href{https://github.com/ipsita-ghosh-1/EDG-IRLS}{github\_SEGRED}

Main paper: 
 introduction: 
 Euclidean Distance Geometry (EDG) problems have applications spanning diverse domains, from comprehending protein structures through molecular conformations [MMLL19, MMWL22] , prediction of molecular conformations in computational chemistry [DPRV15, SWBB97] and aiding dimensional reduction in machine learning [TSL00] to facilitating localization in sensor networks [FQX19] , coupled with its role in solving partial differential equations on manifolds [LL17] . This emphasizes its broad impact in computational sciences. [LLMM14] proposes embedding entities based on Distance Geometry Problems (DGP), where object positions are determined based on a subset of pairwise distances or inner products which significantly reduces computational complexity compared to traditional word embedding methods.Problem 1. Mathematically, consider a collection of n points p i in an r-dimensional Euclidean space with coordinates P = [p 1 , p 2 , ..., p n ] ∈ R r×n whose pairwise squared Euclidean distances are given by d 2 ij = ||p ip j || 2 for each 1 ≤ i ̸ = j ≤ n where || • || denotes the Euclidean norm. Given only partial information of the {d ij } such that only a subset of cardinality m < n(n -1)/2 is known, the goal is to reconstruct the geometry of the points, that is to recover the point coordinates P.If all the distances between the points are provided, the problem is known as multidimensional scaling [LT24, GCKG23] , and closed solution formula exists. However, this is not the case for the incomplete setup which is the focus of this work, and in which only partial information of the pairwise distances is available.For instance, AlphaFold [SEJ + 20, JEP + 21] has shown effectiveness in predicting the threedimensional structure of protein given as input its amino acid sequence. AlphaFold2, which uses an attention mechanism-based transformer architecture [VSP + 17], is trained on known sequences and structures from the Protein Data Bank [HMBFGG + 00] and determines the distances between the C α atoms of all residue pairs in a protein. Subsequently, as a subproblem, this distance information is used to predict the protein's structure by identifying the foldings within the protein molecule. In the context of this subproblem, where the input is the predicted pairwise distances, a linear mapping can be used to predict geometric coordinates. However, many of the predicted distances might not be accurate, which is why low confidence regions of the resulting structures, as indicated by the predicted local-distance difference test (pLDDT) [JEP + 21], could be masked, and a new structure could subsequently be re-computed based on the distance entries of the medium-to-high confidence regions using high-accuracy solution algorithms for Problem 1.In the context of the Problem 1, we have access to a subset of entries of this D = [d 2 ij ] matrix, defined by the index set Ω ⊂ {1, . . . , n} 2 . Now, we formulate the Gram matrix X 0 = P ⊤ P, residing in the set of real-valued symmetric matrices S n . This matrix has a lower rank r compared to the symmetric distance matrix D ∈ S n = {X ∈ R n×n : X = X ⊤ }, which has a rank of r + 2 [DPRV15] .To make the solution translation invariant, the centroid of the points must be the origin, i.e., n ℓ=1 p ℓ = 0. To ensure this, the Gram matrix should satisfy also satisfy X 0 • 1 = 0, where 1 ∈ R n is the vector of all ones. Based on these observations, we can frame Problem 1 as an instance of the computationally challenging NP-hard rank minimization [Rec11] problem defined by min X∈Sn,X•1=0,X⪰0 rank(X) subject to X i,i + X j,j -2X i,j = D i,j for all (i, j) ∈ Ω, (1) incorporating also the positive semi-definiteness constraint X ⪰ 0, which comes from the observation that X 0 = P ⊤ P is also PSD.From an optimization perspective, this rank-minimization problem is highly complex due to its nonconvexity and non-smoothness. Being an NP-hard problem [TL18, Hen95, MW97] , it is extremely difficult to solve directly. Consequently, a significant amount of existing research [Rec11, DR16, CR09, CW18, Gro11, TL18] has focused on minimizing the convex envelope of the rank function, known as the nuclear norm, instead. However, it was observed that for related low-rank matrix completion problem, such convex relaxations are not as data-efficient as other formulations [ALMT14, TW13] , while posing also computational limitations [CC14] .In this paper, we study the question of, given incomplete distance information Problem 1, how many samples are necessary to ensure accurate recovery? [TL18] established convergence with O(νrn log 2 n) uniformly random samples for solving (1) using a nuclear norm surrogate. While there are numerous results for non-convex methods in other low-rank recovery problems [KV21, CLC19, Van13] , there are currently only few non-convex algorithm specifically designed for EDG problems available [TL18, SLCT23, ZCZ22, NKKS19, LS24]. To the best of our knowledge, only [NKKS19, LS24] provide convergence guarantees in the non-convex setting. Unlike generic lowrank matrix recovery problems [RFP10] , EDG problems share the complexities of low-rank matrix completion problems, such as the absence of a direct uniform null space property [RXH11] or a restricted isometry property [RFP10] . Additionally, the underlying measurement basis in EDG problems is not orthonormal, making it impossible to directly use the analysis of standard matrix completion to provide guarantees for techniques using the Gram matrix-based low-rank modelling (1) of Problem 1.In this paper, we propose and analyze an algorithm for the EDG problem based on the iteratively reweighted least squares framework (MatrixIRLS, see Section 3), for which we show that m = Ω(νrn log(n)) (where ν is the coherence factor) randomly sampled distances are sufficient to guarantee local convergence to the ground truth, X 0 with a quadratic rate as shown in Theorem 4.3, from which the geometry of the points P is trivially recovered. The sample complexity assumption of Theorem 4.3 matches the lower bound of low-rank matrix completion problems as established in [CT10] . In Section 4, we construct a dual basis (Lemma 4.4) that spans S n which enables us to show the restricted isometry property (RIP) restricted to a tangent space of the manifold of symmetric rank-r in Theorem 4.5 for our approach, which can be of independent interest for the analysis of other nonconvex algorithms.While the convergence statement of Theorem 4.3 only applies in a local neighborhood of a ground truth, the indicated data-efficiency of the proposed method is numerically validated through different experiments in Section 5 on synthetic and real data in comparison to the state-of-the-art methods. Furthermore, we demonstrate that MatrixIRLS method is robust to ill-conditioned data, further highlighting its flexibility. In the Supplementary material, we discuss the limitations in Appendix B. Additionally, we provide proofs related to the the theoretical results of Section 4 in Appendices B and C. We further discuss the numerical considerations of our experiments in Appendix E. We discuss about the computational complexity of the algorithm in detail in Appendix F. 
related work: 
 Early research on EDG problems focused on establishing its mathematical properties like defining the conditions under which a distance matrix can be represented in Euclidean space [Gow85, YH38, YH38, Alf05] .A comprehensive overview of EDG applications, including molecular conformations, wireless sensor networks, robotics, and manifold learning, is available in [Lib20] . Motivated by the molecular conformation problem, [BJ95] relates this Euclidean distance matrix completion to a graph realization problem by showing that such matrices can be completed if the graph is chordal. Other than graph theoretic approaches, as a technical tool, various optimization strategies [Hen95, MW97, Tro00] have been deployed to solve EDG problems. [AKW99] proposes a primaldual interior point algorithm that solves an equivalent semi-definite programming problem. However, none these works provide theoretical reconstruction guarantees in the incomplete setup of Problem 1.While providing an accurate modelling of Problem 1, the rank minimization formulation (1) poses challenges due to its non-convex and non-smooth nature. There is a mature existing literature [DR16, CR09, CW18, Gro11] around replacing the rank function, rank(X), with the sums of its singular values σ i (X) (also known as the nuclear norm). Building on the rank minimization formulation (1) of the EDG problem, [TL18] minimizes the convex nuclear norm surrogate of the inferred Gram matrix. They propose a dual basis approach that enables a theoretical guarantee for this type of nuclear norm minimization formulation of the EDG problem.From a practical point of view, it is well-known that the convex approach is computationally intensive, as the arithmetic complexity is cubic in n, the dimension of X 0 [CC14] . It also tends to demand more data samples than non-convex alternatives, making it less efficient in terms of data [ALMT14, TW13] .To mitigate these issues, recent studies have shifted focus to non-convex methods such as matrix factorization [SL16, MWCC20, ZCZ22] . These methods optimize a function involving a data-fit objective regularized by squared Frobenius norms of the factor matrices, which are computationally more feasible and data-efficient. Few "non-convex" algorithms are based on matrix factorization like the work in [BM03] .Based on a similar formulation, ScaledSGD [ZCZ22] is a preconditioned stochastic gradient descent method aimed at robustness for ill-conditioned problems. Additionally, some of the most effective techniques for low-rank matrix completion involve minimizing smooth objectives on the Riemannian manifold of fixed-rank matrices, providing scalability and the potential to reconstruct the matrix with fewer samples, although they lack strong performance guarantees [Van13, WCCL20, BA15, BNZ21] . ReiEDG [SLCT23] is a Riemannian-based gradient descent strategy utilizing the sampling operator on Ω, which is non-convex approach to solving the EDG problem. However, it does not provide convergence guarantees. To the best of our knowledge, [NKKS19, LS24] are the non-convex approaches for solving the EDG problem in the Gram-matrix-based low-rank modeling (1) of Problem 1. The work in [NKKS19] proposes an algorithm based on Riemannian optimization over a manifold and provides convergence guarantees. These guarantees are derived from extended Wolfe conditions. However, it is not explicitly detailed how these convergence guarantees depend on problem parameters such as the sampling model and the number of samples (see Remark III.8 in [NKKS19] ). Similarly, the study in [LS24] employs a Riemannian framework and provides local convergence guarantees under Bernoulli sampling. Nonetheless, [LS24] does not clarify whether the proposed algorithm achieves local linear convergence. In contrast to these studies, the algorithm proposed in this paper achieves local quadratic convergence under uniform sampling of the distances.To handle the non-smoothness and non-convexity of rank minimization problems of the type (1), iteratively reweighted least squares (IRLS) algorithms take a different route than methods based on [BM03] or Riemannian methods by minimizing a sequence of quadratic majorizing functions of smoothed rank surrogates [DDFG10, FRW11, MF10, KS18, KV21]. IRLS algorithms have been studied extensively over the years, as indicated by [Law61, GR97, Bur12] . In the context of lowrank matrix completion, IRLS algorithms are known to be among the most data-efficient methods available, while being amenable to a rigorous convergence analysis [FRW11, MF10, KS18, KV21]. Most recently, [KV21] provided an improvement on previous instantiations of the IRLS framework [FRW11, MF10, KS18] for low-rank optimization problems by providing an improved reweighting strategy, for which the authors show a local convergence guarantee that is applicable for low-rank matrix completion, given random entrywise samples of minimal sample complexity. The algorithm we propose in Section 3 is similar to [KV21, Algorithm 1] and Theorem 4.3 follows partially the proof strategy of a related result in [KV21] . However, the setup of Problem 1 does not allow a direct adaptation of both the implementation and analysis of [KV21] due to the non-orthogonality of the measurement basis. 
matrixirls for euclidean distance geometry: 
 In this section, we provide a detailed outline and description of the iteratively reweighted least squares method in the context of the EDG reconstruction problem. To this end, we define preliminaries for stating the algorithm. MatrixIRLS, defined in Algorithm 1 below, can be interpreted as a hybrid of a smoothing method [Che12] and a Majorization-Minimization algorithm [SBP17] . In particular, the proposed algorithm minimizes smoothed log-det objectives defined asEQUATIONwhich is a continuously differentiable function with ϵ -2 -Lipschitz gradient [KV21] .We can decompose X by X = U dg(γσ)U ⊤ , where γ ∈ {+1, -1}, It is clear that the optimization landscape of F ϵ crucially depends on the smoothing parameter ϵ. Instead of minimizing F ϵ directly, our method minimizes, for k ∈ N, ϵ k > 0 and X (k) a quadratic model that is related to the secondorder Taylor expansion of the function f ϵ at the current iterate and its information is encoded in a weight operator [KV21] defined below in Definition 3.1.). X ∈ S n be a matrix with singular value decomposition X (k) = U dg(γσ)U ⊤ , where γ ∈ {+1, -1} i.e., U ∈ S n are orthonormal matrices. We define the weight operator core matrix H X,ε ∈ S n of X for smoothing parameter ε > 0 such that (H X,ε ) ij := max(σ i , ϵ) max(σ j , ϵ) -1 for each i, j ∈ {1, . . . , n} and the weight operator W X,ε :EQUATIONwhere • denotes the entrywise ore Hadamard product of two matrices.Our method is designed to provide iterates that satisfy the constraints of the formulation (1) at each iteration. They can be encoded using the following definition.Definition 3.2. Given n ∈ N, let I = {α = (i, j) | 1 ≤ i < j ≤ n} be the index set of upper triangular indices.We define the operator basis {w α } α∈I∪{(i,i):i∈{1,...,n}} whereEQUATIONwhere e i ∈ R n is the i-th standard basis vector.Given a set (or multiset) of indices Ω = {α 1 , . . . , α m } ⊂ I of cardinality m = |Ω|, we define the measurement operator A = A Ω : S n → R m+n which maps X ∈ S n to A(X) whose ℓ-th coordinate is defined as A(X) ℓ = ⟨w α ℓ , X⟩ for ℓ ≤ m and as A(X) ℓ = ⟨w (ℓ-m,ℓ-m) , X⟩ for ℓ > m.The basis {w α } α of Definition 3.2 can be considered as an extended definition of the primal basis used in [TL18, LT24] , but additionally is able to encode the constraint X • 1 = 0 which guarantees that the Gram matrix corresponds to points whose centroid is located at the origin. Accordingly, we can define the constraint set corresponding to Gram matrices of points that are centered and simultaneously satisfy the pairwise distance constraints of Problem 1 as {X ∈ S n : A(X) = [D Ω ; 0]}. The algorithm below minimizes the quadratic, majorizing model of the objective F ϵ k (•) given k, while satisfying the measurement operator based on the sampled distances. Equivalently, we can define the main computational step of the method asEQUATIONwhere W k : S n → S n is defined as W k := W X (k) ,EQUATIONUpdate smoothing: Compute r + 1-th singular value of X (k) to updateEQUATIONUpdate weight operator: Forr k := |{i ∈ [n] : σ i (X (k) ) > ϵ k }|, compute the first r k singular values σ (k) i := σ i (X (k)) and matrices U (k) ∈ R n×r k with leading r k left singular vectors of X (k) to update W k = W X (k) ,ϵ k defined in Definition 3.1. end forThe computational complexity of this above algorithm can be computed from the steps (6) and (7). In our numerical implementation, we largely follow the tangent space formulation of the weighted least squares step (6), cf. [KV21, Section 3], which involves the solution of an order O(nr k ) = O(nr) linear system if r = r is chosen as the ground truth rank. An additional difficulty we overcame in the provided reference implementation arises from the fact that AA * is not the identity. The per-iteration time complexity of our method is dominated by O((mr + r 2 n)N 0 inner ), where N 0 inner is the number of inner iteration bound of the iterative linear system solver. Detailed FLOPs calculation is shown in Appendix F. matrix is concentrated in few entries, and this observation is typically captured by the notion of incoherence [CR09, Rec11, Che15] . For the purpose of the EDG problem of interest, we use the following coherence notion, which has appeared in similar form in [TL18] .Definition 4.1 (Coherence for Gram matrices in the EDG problem, [TL18] ). Let X ∈ S n be of rank r. Let T = T X = {ZX + XZ ⊤ : Z ∈ R n×n } be the tangent space onto the rank-r manifold M r = {Z ∈ S n : rank(Z) = r} at X. We say that X has coherence ν with respect to the basis {w α } α∈I of the subspace {X ∈ S n :X1 = 0} if max α∈I β∈I ⟨P T w α , w β ⟩ 2 ≤ 2ν r n and max α∈I β∈I ⟨P T v α , w β ⟩ 2 ≤ 4ν r n ,where P T : S n → S n denotes the projection operator onto T and {v} α∈I is a dual basis of {w α } α∈I , which means that ⟨w α , v β ⟩ = δ α,β for each α, β ∈ I ( δ α,β = 1 for α = β and equal to 0 otherwise).Remark 4.2. In [TL18, Definition 1], the coherence constant ν was required to satisfy a third condition (see [TL18, (Ineq. 14)]). However, this condition is not needed for our proofs, which is why we can use the weaker definition of Definition 4.1. Similar improvements for the standard basis incoherence notion were achieved in [Che15] . In [TL18, Lemma 21], it was shown that up constants, the definition above is equivalent to a coherence condition with respect to the standard basis [Rec11, Che15] .Following the conventional sampling approach in the existing literature [Rec11, Che15, CWB08], the index set is Ω = (i ℓ , j ℓ ) m ℓ=1 ⊂ I contains m samples drawn uniformly at random without replacement. Theorem 4.3 (Local convergence of MatrixIRLS for EDG with Quadratic Rate). Let X 0 ∈ S n be a matrix of rank r that is ν-incoherent, and let A : S n → R m+n be the measurement operator corresponding to an index set Ω ⊂ I of size m = |Ω| that is drawn uniformly without replacement. There exist constants C * , C and C such that the following holds. (a) If the sample complexity fulfills m ≥ Cνrn log n, and if (b) the output matrix X (k) ∈ S n of the k-th iteration of MatrixIRLS for EDG with inputs y = A(X 0 ) and r = r updates the smoothing parameter in (7) such thatϵ k = σ r+1 (X (k) ) and fulfills ∥X (k) -X 0 ∥ S∞ ≲ m 3 2 C * κL 2 (log n) 3 2 √ n-r σ r (X 0 ) where κ = σ 1 (X 0 )/σ r (X 0 )is the condition number of X 0 , then the local convergence rate is quadratic in the sense that ∥X(k+1) -X 0 ∥ S∞ ≤ min(µ∥X (k) -X 0 ∥ 2 S∞ , ∥X (k) -X 0 ∥ S∞ ) with µ = m C 2 L log n 1 4(1+6κ)σr(X 0 )and furthermore X (k+ℓ) ℓ→∞ ---→ X 0 with high probability.(The values of the constants C * = 10 5 , C = 21 √ 5,C = 4900 are explicitly derived in the Supplementary material.)In other words, Theorem 4.3 indicates Algorithm 1 converges to the ground truth with high probability with a sample complexity of Ω(νrn log n), if initialized close to the ground truth Gram matrix. We refer to Appendix D for its proof. Our theorem's sample complexity requirement aligns with the lower bound for generic low-rank matrix completion problems, as established in [CT10] . We note that this theorem only provides a local convergence guarantee for Algorithm 1. This is in line with the strongest known results for IRLS algorithms optimizing non-convex objectives [For10, KV21, PKV22]. We provide numerical evidence in Section 5 that the minimal sample complexity assumption (a) of Theorem 4.3 indeed captures the generic reconstruction ability of the method.To the best of our knowledge, Theorem 4.3 represents the first convergence guarantee for any algorithm for Problem 1 that applies at the optimal order Ω(νrn log n) of provided pairwise distances, and furthermore, the first theoretical guarantee for any non-convex optimization framework for Problem 1. 
dual basis construction and local restricted isometry property on tangent spaces: 
 While a convergence result similar to Theorem 4.3 had been previously obtained for an IRLS algorithm for low-rank matrix completion [KV21] , an adaptation of the proof of [KV21] to the EDG setting is not possible due to non-orthogonality of the basis {w α } α of Definition 3.2.In order to prove Theorem 4.3, we establish a restricted isometry property of a suitably defined sampling operator (see (9)) with respect to the tangent space of the manifold of symmetric rank-r matrices at the ground truth Gram matrix X 0 = P ⊤ P. To formulate this sampling operator and the respective RIP condition, we construct a dual basis to the measurement basis {w α } α of Definition 3.2. Lemma 4.4 (Dual Basis Construction). Let n ∈ N, I = {(i, j) | 1 ≤ i < j ≤ n} be the index set and {w α } α∈I∪I D be the primal basis of Definition 3.2 withEQUATIONwhere a i = e i -1 n 1 for i ∈ {1, . . . , n}, then {v α } α∈I∪I D is a dual basis with respect to {w α } α∈I∪I D , i.e., {v α } α and {w α } α are bi-orthogonal.Lemma 4.4 extends the dual basis construction of [TL18, LT24] , in which the duality of {v α } α∈I with respect to {w α } α∈I was shown. The proof of Lemma 4.4 is detailed in Appendix B.1.Unlike the basis pair of [TL18, LT24] , our bases span the entire space of symmetric matrices S n (see Appendix B.2), which enables us to show the following restricted isometry property.Q Ω ). Let L = n(n -1)/2, 0 < ϵ ≤ 1 2 ,and Ω ⊂ I be a multiset of size m sampled independently with replacement. Define the sampling operatorQ Ω : S n → S n such that Q Ω (X) := L m α∈Ω∪(i,i) n i=1 ⟨X, w α ⟩v α (9)where w α and v α as in (4) and (8), respectively. Let X 0 ∈ S n be a ν-incoherent matrix whose tangent space onto the manifold M r of symmetric rank-r matrices is denoted as T 0 = T X 0 . Let P T0 : S n → S n be the projection operator associated to T 0 . Then P T0 Q * Ω P T0 -P T0 S∞ ≤ ϵ holds with probability at least 1 -2 n provided that m ≥ (49/ϵ 2 )νnr log n.The dual basis as discussed in Section 3 along with the extension for the diagonal entries together spans the space of n × n symmetric matrices. This construction has been crucial in proving the RIP restricted to the tangent space T X of rank constrained smooth manifold M r . This construction could also be valuable for analyzing other non-convex algorithms. A detailed proof of Theorem 4.5 is provided in Appendix C This proof is achieved by using concentration inequalities like the Matrix Bernstein inequality theorem C.1 in multiple lemmas stated and proved in detail in the supplemental material. Since the existing literature for nonconvex approaches for solving Problem 1 lacks this property, this approach of establishing RIP by restricting it to the tangent space can be useful for the analysis of other nonconvex methods. The RIP condition, originally introduced in the context of compressed sensing [Can08] , is a fundamental assumption in the literature on low-rank matrix recovery ([ZL12, Che15, KV21, CZ13]). A tangent-space restricted RIP has been has been useful for analyzing the convergence and performance properties of other non-convex methods [GLM16, LZ23, LHLZ20, ZLTW18] in a non-EDG setting.We count the degrees of freedom in the spectral decomposition of a rank-r symmetric matrix: A = UDU ⊤ where U is an orthogonal matrix and D is a diagonal matrix consisting of the eigenvalues of A on the diagonal. In U there are total r unit norm constraints and there are total r(r-1) 2constraints to the orthogonality of the columns vectors of U which follow from ⟨u i , u j ⟩ = 0 for i ̸ = j and i ∈ {1, 2, ..., r} and j ∈ {1, 2, ..., r}. Therefore, the total number of constraints isr(r-1) 2 + r = r(r+1) 
numerical experiments: 
 We evaluate the performance of MatrixIRLS, Algorithm 1, for instances of the EDG reconstruction Problem 1 in terms of data efficiency across multiple datasets in comparison to other state-of-the-art methods in the literature. We compare the performance of MatrixIRLS with three other algorithms: (a) ALM, an augmented Lagrangian method that minimizes the non-convex formulation defined by min P∈R r×n Tr(P ⊤ P) subject to R Ω (P ⊤ P) = R Ω (X 0 ), with R Ω being equal to the sampling operator Q Ω restricted to Ω, which is based on a Burer-Monteiro factorization of the Gram matrix, and which has been studied in the numerical experiments of [TL18] , (b) ScaledSGD [ZCZ22] which is a preconditioned stochastic gradient descent method designed to be robust with respect to ill-conditioned problems, (c) RieEDG [SLCT23] which is a Riemannian-based gradient descent approach based on the sampling operator (9) restricted to Ω. The choice of these algorithms is based on their robustness to noise as claimed in the respective papers.We first consider a synthetic data, where we select n = 500 points P 0 = [p 1 , . . . , p n ] ∈ R r×n from a standard Gaussian distribution at random such that (p i ) j i.i.d.∼ N (0, 1) for all i, j, which defines the ground truth Gram matrix X 0 = P 0⊤ P 0 . We are provided with m = |Ω| Euclidean distances, where the point index pair set Ω ⊂ {(i, j) ∈ [n] × [n], i < j} is sampled uniformly at random. This is parametrized by the oversampling factor ρ = m nr-r(r-1)/2 where the denominator is the degrees of freedom (discussed in Appendix E.4). To understand the efficiency of the algorithm over a range of ranks r and across different oversampling factors ρ, we conduct phase transition experiments for all the above mentioned algorithms. We define a successful recovery as the case of the relative Procrustes distance d Procrustes (P rec , P 0 ) (discussed in Appendix E.3) between the recovered matrix P rec and ground truth coordinate matrix P 0 does not exceed a tolerance threshold tol rec = 10 -3 . We chose the Procrustes distance at it is a shape preserving distance that accounts also for differences in scaling and alignment between the reconstructed geometries. We observe the performance of the algorithms as shown in the Figure 1 for ranks between 2 to 5 and a oversampling factor ranging from 1 to 4 over 24 instances.In Figure 1 , each entry on the figure represents the probability of success of an algorithm for the given rank ground truth rank r and oversampling factor ρ over 24 instances. 1 In terms of the recovery of the ground truth, we notice a comparable performance for the MatrixIRLS and ALM. However, the other two algorithms RieEDG and ScaledSGD, for the given success tolerance, the recovery of a ground truth is only possible if more samplies corresponding to a larger oversampling factor are provided, for any of the ranks r considered. This emphasizes that MatrixIRLS is able to achieve state-of-the-art data efficiency. In order to understand the scalability of the proposed algorithm, we have provided the time of completion of MatrixIRLS in Table 1 . The first two rows shows that for fixed oversampling factorρ = 3 , we are able to recover the points in the magnitude of 10 4 , is just 13.7 minutes. So to further test the scalability, we also looked at the time to completion when less number of samples are provided ( ρ = 2.5). In that setup n = 10000 takes 57.5 minutes to recover with high precision.We now assess the performance of the different EDG methods on ill-conditioned data. Ill-conditioning in the point matrix in particular arises in situations where, for example, r-dimensional points are approximately following an r -1-dimensional geometry, a situation which often arises when the available distance information if affected by outliers [ZCZ22] .Similar to the setup above, we construct Gram matrices X 0 ∈ S n of rank r = 5 with condition number κ = σ 1 (X 0 )/σ r (X 0 ) = 10 5 corresponding to n = 400 random data points P 0 ∈ R r×n generated with random orthogonal singular vectors and singular values σ i (X 0 ) interpolated between κ and 1 with decay of order O(i -2 ). Figure 2 shows the success probability of different algorithms for ill-conditioned ground truths.To have a deeper insight into the algorithm's performance, Figure 3 shows the per-iterate relative reconstructed error of the four algorithms on this ill-conditioned data at an oversampling factor of ρ = 2 for a representative problem instance. This reconstructed error refers to the Procrustes distance between the ground truth P 0 and the recovered P k at each iteration k. We observe that for ill-conditioned data, MatrixIRLS is the only method that can recover the ground truth up to a reasonable precision, achieving a relative error of ≈ 10 -8 after around 35 iterations (top of Figure 3 ), whereas the other methods do not achieve errors below 10 -3 even after 100000 iterations. While one iteration of MatrixIRLS typically has a longer runtime than an iteration of any of the other algorithms due to its second-order nature, we also provide a visualization of the observed runtimes of the methods in the bottom of Figure 3 . It can be seen that MatrixIRLS is able reconstruct the geometry of the challenging problem in around 200 seconds up to a high precision. Additionally, we run a study on the algorithms' behavior across different oversampling factors between ρ = 1 and ρ = 4 for 24 random problem instances, visualized in Figure 4 . The box plots indicate visualize median relative Procrustes error together with the relevant 25% and 75% quantiles of the observed error distribution. We observe that MatrixIRLS has consistent convergence for ill-conditioned data even for lower oversampling rates as long as ρ ≥ 1.5. To assess the performance of Algorithm 1 in realistic setups, we consider the task of molecular conformation, i.e., we aim to reconstruct the 3D structure of a molecule from partial information about the distances of its atoms. For our experiments, we determine the structures of a protein molecule (1BPM) from the Protein Data Bank [HMBFGG + 00], which is collected using X-ray diffraction experiments or nuclear magnetic resonance (NMR) spectroscopy.The goal of this experiment is to reconstruct the point matrix P 0 from |Ω| distance samples as defined in Problem 1. For the 1BPM protein data the rank of the input matrix is 3. We conduct the experiments corresponding to oversampling factors ρ between 1 and 4. Like the previous setup, here successful recovery refers to the case where the relative Procrustes distance, that is the metric d Procrustes (P rec , P 0 ) between the recovered matrix P rec and ground truth coordinate matrix P 0 does not exceed a tolerance threshold tol rec across 24 independent realizations. The error analysis for this experiment is shown in Figures 7a and 7b in Appendix E.It is evident that MatrixIRLS and ALM have comparable results in recovering the geometry of the data given lesser samples, where as algorithms like ScaledSGD or RieEDG are unable to reconstruct geometries even when the oversampling rate is as high as ρ = 4. In the fig. 5 , we see a convergence with high precision for MatrixIRLS from ρ ∼ 2.5 for Protein which means that it recovered the ground truth with around 0.5% samples. In the Figure 6 the reconstructed protein in blue is aligned with the ground truth structure of the protein in pink. We can see that when the oversampling is 3.5, which is when there is 0.6% samples available, the reconstruction exactly matches with the ground truth.In order to understand the runtime of the different algorithms Table 2 , reports the reconstruction times for each of the four algorithms applied to the 1BPM Protein data (with datapoints n = 3672) in a low-data regime with oversampling factor of ρ = 3. . It can be seen that MatrixIRLS is with 7.08 minutes around 3 times faster than ALM, which needed 23.04 minutes until convergence. While these times certainly depend on the precise choice of stopping criteria for each algorithm (we used the ones indicated in Appendix E ), this shows that the proposed method is competitive in terms of clock time. For the implementation of the other algorithms, we use the authors' code for the respective approaches (discussed in Appendix E). We include another set of experiments on US cities data [UU20], in the Appendix E. 
conclusion: 
 In this paper, we address the challenge of reconstructing suitable geometric configurations using minimal Euclidean distance samples. By leveraging continuous and non-convex rank minimization formulations, we develop a variant of the iteratively reweighted least squares (IRLS) algorithm and establish a local convergence guarantee under the condition of a minimal random set of observed distances. Our contribution also includes the proof of a restricted isometry property (RIP) restricted to the tangent space of the manifold of symmetric rank-r matrices, a result which might be of independent interest for the analysis of other non-convex methods. As future work further analysis on the global convergence can be established. Through numerical validation we conclude that the algorithm is able to achieve state-of-the-art data efficiency. 
broader impacts: 
 Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?Answer: [Yes] Justification: In the introduction and in related work we highlight the impact of the work. These provide foundational algorithmic research to reconstruction algorithms which has applications in molecular conformation, sensor network localization. There is no potential negative docial impact of this work.• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations. • The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?Answer: [NA]Justification: The paper uses data that are publicly available. The data sources are cited properly in Section 5.Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. • We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort.12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?Answer:[Yes]Justification: Throughout the paper we cite the different information that we used from the existing literature. Additionally, in the Appendix E, we provide links to the Github repositories that we used for reference.Guidelines:• The answer NA means that the paper does not use existing assets.• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a URL. • The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided. • If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset. • For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. • If this information is not available online, the authors are encouraged to reach out to the asset's creators. 13. New Assets Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets? Answer: [Yes] Justification: The link to the anonymized version of the repository containing our experiments have been provided in Appendix E. Guidelines:• The answer NA means that the paper does not release new assets.• Researchers should communicate the details of the dataset/code/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc. • The paper should discuss whether and how consent was obtained from people whose asset is used. • At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file. 
crowdsourcing and research with human subjects: 
 Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)? Answer: NA Justification: The paper does not involve crowdsourcing nor research with human subjects.r Guidelines:• The answer NA means that the paper does not involve crowdsourcing nor research with human subjects. • Including this information in the supplemental material is fine, but if the main contribution of the paper involves human subjects, then as much detail as possible should be included in the main paper. Guidelines:• The answer NA means that the paper does not involve crowdsourcing nor research with human subjects. • Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper. • We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution. • For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review. 
 
 Appendix: For our method, the limited convergence radius leads to the convergence guarantee locally. Additionally, the runtime per iteration is larger compared to other non-convex methods [ZCZ22, TL18, SLCT23] but our method achieves convergence in less than a tenth of the number of iterations of the other methods. For some datasets like the Protein (1BPM) [HMBFGG + 00], ALM has a superior performance for oversampling factor (ρ) between 1 and 1.5.In this section, we provide the proof the bi-orthogonality of the basis defined in Lemma 4.4 with respect to the measurement basis Definition 3.2, and furthermore establish properties of the basis pair in Appendix B.2 that will be useful later on.Proof of Lemma 4.4. For simplicity of the proof, we denote the two cases of Definition 3.2, separately. This helps us to divide the proof into multiple subcases.With L = n(n-1) 2 , let V = v (1,2) , . . . , v (n-1,n) ∈ R n 2 ×Lbe a matrix whose columns enumerate the vectorized dual basis elements v α where α = (i, j) ∈ I, and letV E ∈ R n 2 ×n consist of vectorized dual basis vectors v (i,i) with i ∈ {1, . . . , n}. Similarly, W ∈ R n 2 ×L andW E ∈ R n 2 ×n are defined using vectorized basis elements {w α } α∈I∪I D . As establishing duality between the two basis sets is equivalent to bi-orthogonality of the basis elements, it remains to show that the extended basis matricesṼ = [V|V E ] ∈ R n 2 ×L+n and W = [W|W E ] ∈ R n 2 ×L+n satisfy the relationship W⊤ Ṽ = W ⊤ V W ⊤ V E W ⊤ E V W ⊤ E V E = Id,where Id is the identity matrix.In coordinate-wise notation, this means that for α = (i, j) and β = (k, l), with 1 ≤ i ≤ j ≤ n and 1 ≤ k ≤ l ≤ n, we must show that ⟨ wα , ṽβ ⟩ = δ α,β where wα and ṽβ are columns of W and Ṽ respectively.From Claim 3.1 of [LT24] , it is clear that V is the dual of W. We require to show the biorthogonality for the extended parts.First we show that, for α of the form (i, i),EQUATIONLet us decompose the summand ⟨e i • 1 ⊤ , a i a ⊤ i ⟩ in detail, so that we can derive the remaining calculations in a similar manner,⟨e i • 1 ⊤ , a i a ⊤ i ⟩ = ⟨e i • 1 ⊤ , e i e ⊤ i ⟩ -⟨e i • 1 ⊤ , 1 n 1e ⊤ i ⟩ -⟨e i • 1 ⊤ , 1 n e i 1 ⊤ ⟩ + ⟨e i • 1 ⊤ , 1 n 2 1 • 1 ⊤ ⟩ = 1 - 1 n -1 + 1 n = 0.Putting this value back in (10), we can complete the rest of the calculation in a similar manner and finally arrive at the following,⟨w α , v α ⟩ = ⟨w i,i , v i,i ⟩ = 1 2 (⟨e i • 1 ⊤ , e i e ⊤ i ⟩ + ⟨1 • e ⊤ i , e i e ⊤ i ⟩ -⟨e i • 1 ⊤ , a i a ⊤ i ⟩ + ⟨1 • e ⊤ i , a i a ⊤ i ⟩ = 1 2 (1 + 1 -0 -0)= 1. Now, let us look at the case where α = (i, i), and β = (j, j) and i ̸ = j. Using a similar calculation we can show that,⟨w α , v β ⟩ = ⟨w i,i , v j,j ⟩ = 1 2 ⟨e i • 1 ⊤ -1 • e ⊤ i , e j e ⊤ j -a j a ⊤ j ⟩ = 1 2 (⟨e i • 1 ⊤ , e j e ⊤ j ⟩ -⟨1 • e ⊤ i , e j e ⊤ j ⟩ -⟨e i • 1 ⊤ , a j a ⊤ j ⟩ + ⟨1 • e ⊤ i , a j a ⊤ j ⟩ = 1 2 (0 -0 -0 + 0)= 0. Also, between the w α and v β where α = (i, j) with 1 ≤ i < j ≤ n and β = (k, k) with 1 ≤ k ≤ n, we can show a similar relation by⟨w α , v β ⟩ = ⟨w i,j , v k,k ⟩ = 1 2 ⟨e i,i + e j,j -e i,j -e j,i , e k e ⊤ k -a k e ⊤ k ⟩ = 1 2 δ i,k + δ j,k -0 -0 -(δ i,k - 1 n )(δ i,k - 1 n ) -(δ j,k - 1 n )(δ j,k - 1 n ) + 1 2 (δ i,k - 1 n )(δ j,k - 1 n ) = 1 2 δ i,k + δ j,k -δ i,k δ i,k + 1 n δ i,k + 1 n δ j,k - 1 n 2 -δ j,k δ j,k + 1 n δ j,k + 1 n δ j,k + 1 2 - 1 n 2 + 2δ i,k δ j,k -2 1 n δ j,k -2δ i,k + 2 n 2 = 1 2 (2δ i,k δ j,k ) = δ i,j = 0.Finally, we need to show for w α with α = (i, i) for 1 ≤ i ≤ n and v β where β = (k, l) with 1 ≤ k < l ≤ n that the respective dot product is 0. This is verified by the calculation⟨w α , v β ⟩ = - 1 4 ⟨e i • 1 ⊤ + 1 • e ⊤ i , e k e ⊤ l + e l e ⊤ k ⟩ = - 1 4 (⟨e i • 1 ⊤ , e k e ⊤ l ⟩ + ⟨e i • 1 ⊤ , e l e ⊤ k ⟩ + ⟨1 • e ⊤ i , e k e ⊤ l ⟩ + ⟨1 • e ⊤ i , e l e ⊤ k ⟩) = - 1 4 (δ ik δ il + δ ik δ il + δ ik δ il + δ ik δ il ) = - 1 4 (δ ik δ il + δ ik δ il + δ ik δ il + δ ik δ il ) = - 1 4 (4δ ik δ il ) = -δ ik δ il = -δ kl = 0Hence combining the above parts we have proved that Ṽ is the dual of W.Lemma B.1. W Ṽ⊤ maps n × n symmetric matrices to itself.Proof. Let x be a vectorized representation of a symmetric matrix X. Then, using the dual basis expansion, x can be represented as: x = W ṼT x. We now apply the the operator W ṼT to x.( W ṼT ) W ṼT x = W ṼT x.where the last equality follows from the fact that Ṽ is dual to W.Lemma B.2. For orthonormal basis {w α } and {v α } as defined in Definition 3.2 , the spectral norm of W Ṽ⊤ is 1, that isEQUATION)Proof. By definition, ∥ W ṼT ∥ S∞ = max ∥x∥2=1 ∥ W ṼT x∥ 2 . Note that ∥ W ṼT x∥ 2 2 = x T Ṽ( WT W) ṼT x.By construction of the dual basis, ( WT W) = ( ṼT Ṽ) -1 . Therefore, Ṽ( WT W) ṼT is an orthogonal projection operator onto the column space of Ṽ. If x is a vectorized representation of any symmetric matrix, using Lemma B.1 the projection operator is an identity operator. Hence, its operator norm is 1.In this section, we first state the matrix Bernstein [Tro11] concentration inequality for rectangular matrices, which is repeatedly used in the proofs further in the section. We then provide a proof for the bound of our sampling operator Q Ω and then we conclude this section by proving the RIP restricted to the tangent space of manifold of symmetric rank-r matrices at the ground truth X 0 .Theorem C.1 (Matrix Bernstein for rectangular matrices [Tro11, Theorem 1.6]). Consider the finite sequence {Z k } of independent, random matrices with dimension d 1 × d 2 . Assume that each random matrix satisfiesE[Z k ] = 0 and ∥Z k ∥ ≤ R almost surely, Define σ 2 := max k E[Z k Z * k ] , k E[Z * k Z k ]Then, for all t ≥ 0,P k Z k ] ≥ t ≤ (d 1 + d 2 ) • exp -t 2 /2 σ 2 + Rt/3 (12)Lemma C.2 (Spectral Norm Bound for Sampling Operator Q Ω ). Let n be the dimension of the input matrix and let Ω = (i ℓ , j ℓ ) m ℓ=1 ⊂ I be a multiset of double indices fulfilling m < L = n(n -1)/2 that are sampled independently with replacement. Consequently, we have that with probability of at least 1 -2 n 2 , the sampling operator Q Ω : S n → S n of (9) fulfills∥ Q Ω ∥ S∞ ≤ 20L log n m + 1.Proof. We define,R Ω ′ = VS Ω W ⊤ = α∈Ω v α w ⊤ α (13)So from (13),R Ω = L m R Ω ′ Q Ω = R Ω +V E W ⊤ E Q Ω ′ = R Ω ′ + L m V E W ⊤ EIn order to provide a bound for Q Ω , we first evaluate a bound for R Ω . We would use concentration inequality of Theorem C.1, for computing the bound of Q Ω Let us defineEQUATIONWe compute the expectation of z k by,EQUATIONNow let us compute the bound for z kEQUATIONWe use triangle inequality in the first inequality and since the spectral norm of VW ⊤ is 1, we arrive at the second inequality. Now, we can further bound the spectral norm of v α k w ⊤ α k by,||v α k w ⊤ α k || S∞ = max ⟨V, v α k w ⊤ α k W⟩ ||V|| 2 ||W|| 2 ≤ ||v α k || 2 2 ||v α k || 2 ||w α k || 2 2 ||w α k || 2 ≤ 1 • 4 = 4 (17)So from (16), we get||z k || S∞ ≤ 4L + 1 m (18)We take the product of z k and z * k as,z * k z k = ( L m v α k w ⊤ α k - 1 m VW ⊤ ) ⊤ ( L m v α k w ⊤ α k - 1 m VW ⊤ ) = ( L m w α k v ⊤ α k - 1 m WV ⊤ )( L m v α k w ⊤ α k - 1 m VW ⊤ ) = L 2 m 2 (w α k v ⊤ α k v α k w ⊤ α k ) - L m 2 (w α k v ⊤ α k VW ⊤ ) - L m 2 (WV ⊤ v α k w ⊤ α k ) + 1 m 2 (WV ⊤ VW ⊤ )Now taking expectation on both side of the above relation we get;||E(z * k z k )|| S∞ = E L 2 m 2 (w α k v ⊤ α k v α k w ⊤ α k ) + E L m 2 (w α k v ⊤ α k VW ⊤ ) + E L m 2 (WV ⊤ v α k w ⊤ α k ) + E 1 m 2 (WV ⊤ VW ⊤ ) S∞ ≤ E L 2 m 2 (w α k v ⊤ α k v α k w ⊤ α k ) S∞ + E L m 2 (w α k v ⊤ α k VW ⊤ ) S∞ + E L m 2 (WV ⊤ v α k w ⊤ α k ) S∞ + E 1 m 2 (WV ⊤ VW ⊤ ) S∞ = L 2 m 2 1 L L k=1 (w α k v ⊤ α k v α k w ⊤ α k ) S∞ + 1 L L k=1 L m 2 (w α k v ⊤ α k VW ⊤ ) S∞ + 1 L L m 2 (w α k v ⊤ α k VW ⊤ ) S∞ + 1 m 2 ≤ L 2 m 2 max ||w α k v ⊤ α k v α k w ⊤ α k || S∞ + L m 2 WV ⊤ VW ⊤ S∞ + L m 2 WV ⊤ VW ⊤ S∞ + 1 m 2 ≤ 16L 2 + 2L + 1 m 2 (19)Since expectation is a linear operator we can preserve the equality in the first equality, then we use triangle inequality in the first inequality. We further use that the spectral norm of WV ⊤ is 1 and arrive at the second inequality. Further using (17) we bound the expression in the last inequality.As per Theorem C.1 we now compute the product of z k and z * kz k z * k = ( L m v α k w ⊤ α k - 1 m VW ⊤ )( L m v α k w ⊤ α k - 1 m VW ⊤ ) ⊤ = ( L m v α k w ⊤ α k - 1 m VW ⊤ )( L m w α k v ⊤ α k - 1 m WV ⊤ ) = L 2 m 2 (v α k w ⊤ α k w α k v ⊤ α k ) - L m 2 (v α k w ⊤ α k WV ⊤ ) - L m 2 (VW ⊤ w α k v ⊤ α k ) + 1 m 2 (VW ⊤ WV ⊤ )Now taking expectation on both side of the above relation we get;EQUATIONSo, now by Matrix Bernstein Inequality as stated Theorem C.1,E(z k ) = 0 The value of the bound of z k is denoted by R in Theorem C.1. So R in this case is ||z k || S∞ ≤ R = 4L + 1 m Further σ is defined by max[ k E(z k z * k ), k E(z * k z k )]as per Theorem C.1. We have,EQUATIONRecalling the definition of z k from ( 14)EQUATIONwith a probability 1 -2n exp( -t 2 2 σ 2 + Rt). We have arrived at the above inequality from the definition of Q Ω and that the spectral norm of the extended basis of the dual pair,V E W ⊤ E is 1. Since L = n(n-1) 2 so, L 2 ≤ n 2 L2 . We can further simplify the denominator by the following,EQUATIONWe arrive at the first inequality by using 1 m ≤ L m and 4L+1 3m ≤ 6L 3m So using this relation, let us assume that t = 20L log n m . Further let us simply Equation (23), by using this value of t. Since by Theorem 4.3, m ≥ Kn log n, where K = Cνr, so log n m ≤ 1 Kn , hence we can say simplify Equation (23), further by the following,EQUATIONWe arrive at the third inequality bu using the expression for L and using logn m ≤ 1 Kn . Then, we can deduce the following from Equation (24),exp -t 2 2 σ 2 + Rt 3 = exp -20 2 L 2 log n 2m 60L 2 m ≤ exp - 400 120 log n = exp(3.33 log( 1 n )) ≤ exp(log( 1 n 3 )) = 1 n 3So the probability in (21) becomes 2n • 1 n 3 If we use (22), then we can conclude that that the spectral norm of Q Ω is bounded by 20L log n m + 1 . Hence, we can say that the spectral norm of Q Ω is bounded by 20L log n m + 1 with a probability of 1 -2 n 2 . Now we prove the restricted isometry property of the sampling operator Q Ω with respect to the tangent space of the rank-r matrix manifold at the ground truth, as stated in Theorem 4.5.Proof of Theorem 4.5.EQUATIONNow we calculate the spectral norm of z kEQUATIONThe third inequality follows from the coherence bounds as in Definition 4.1 .z * k z k = ( L m P T0 w α k v ⊤ α k P T0 - 1 m P T0 WV ⊤ P T0 ) ⊤ ( L m P T0 w α k v ⊤ α k P T0 - 1 m P T0 WV ⊤ P T0 ) = ( L 2 m 2 P T0 v α k w ⊤ α k P T0 w α k v ⊤ α k P T0 ) -( L m 2 P T0 v α k w ⊤ α k P T0 P T0 WV ⊤ P T0 ) -( L m 2 P T0 VW ⊤ P T0 P T0 w α k v ⊤ α k P T0 ) +1 m 2 P T0 VW ⊤ P T0 P T0 WV ⊤ P T0 We use that P 2 T0 = P T0 in the second equality to further calculate the expectation of the above equality,E(z * k z k ) = E( L 2 m 2 P T0 v α k w ⊤ α k P T0 w α k v ⊤ α k P T0 ) -E(( L m 2 P T0 v α k w ⊤ α k P T0 P T0 WV ⊤ P T0 ) -E(( L m 2 P T0 VW ⊤ P T0 P T0 w α k v ⊤ α k P T0 )) + E( 1 m 2 P T0 VW ⊤ P T0 P T0 WV ⊤ P T0 ) = E( L 2 m 2 P T0 v α k w ⊤ α k P T0 w α k v ⊤ α k P T0 ) - 1 m 2 P T0 VW ⊤ P T0 WV ⊤ P T0 - 1 m 2 P T0 VW ⊤ P T0 P T0 WV ⊤ P T0 + 1 m 2 P T0 VW ⊤ P T0 P T0 WV ⊤ P T0 = E( L 2 m 2 P T0 v α k w ⊤ α k P T0 w α k v ⊤ α k P T0 ) - 1 m 2 P T0 VW ⊤ P T0 P T0 WV ⊤ P T0since W Ṽ⊤ = Id, we can adjust summand at the seocond equality. Defining the random operator z k := L m P T0 w α k v ⊤ α k P T0 , we see that z k from above satisfies z k = z k -1 m P T0 WV ⊤ P T0 and thereforeE[ z k ] = 1 m P T0 WV ⊤ P T0With this definition, and for the fact that positive semi-definite matrices satisfies ∥A -B∥ S∞ ≤ max(∥A∥ S∞ , ∥B∥ S∞ ), we can bound the spectral norm of the expectation in the following wayEQUATIONFor any M ∈ R n×n , it follows from the definition of z kz k (M) = L m ⟨v α k , P T0 (M)⟩P T0 (w α k ) = L m ⟨P T0 (v α k ), M⟩P T0 (w α k ), z * k (M) = L m ⟨w α k , P T0 (M)⟩P T0 (v α k ) = L m ⟨P T0 (w α k ), M⟩P T0 (v α k )and thusz * k z k (M) = L m ⟨P T0 (w α k ), z k (M)⟩P T0 (v α k ) = L 2 m 2 ⟨P T0 (w α k ), ⟨P T0 (v α k ), M⟩P T0 (w α k )⟩P T0 (v α k ) = L 2 m 2 ∥P T0 (w α k )∥ 2 F ⟨P T0 (v α k ), M⟩P T0 (v α k )Now using the incoherence condition such as Definition 4.1 and [TL18, Lemma 22]. we can argue thatE z * k z k S∞ ≤ L 2 m 2 L max ℓ=1 ∥P T0 (w ℓ )∥ 2 F E P T0 v α k v ⊤ α k P T0 S∞ = L 2 m 2 L max ℓ=1 ∥P T0 (w ℓ )∥ 2 F 1 L L ℓ=1 P T0 v ℓ v ℓ ⊤ P T0 S∞ ≤ L 2 m 2 L max ℓ=1 ∥P T0 (w ℓ )∥ 2 F L max ℓ=1 ∥P T0 (v α k )∥ 2 s∞ ≤ L 2 m 2 (2ν r n ) 1 L (8ν r n ) = L m 2 16 ν 2 r 2 n 2 .Now we get from ( 27)EQUATIONSo by Triangle Inequality,EQUATIONNow we need to compute similar bounds forEQUATIONSo similar to (29) we apply Triangle Inequality on (30) and get,EQUATIONNow, we can approximate this bound in the following way, Since L = n(n-1), so L ≤ n 2 2 this bound can be written asm k=1 E(z k z * k ) S∞ ≤ L m 16 ν 2 r 2 n 2 ≤ 8 ν 2 r 2 m (32)Comparing the equations ( 31) and (32), we can argue that the bound on (31) is larger than that of (32) as L 2 is larger than νr. SoEQUATIONSo, now by Matrix Bernstein Inequality as stated Theorem C.1,E(z k ) = 0 The value of the bound of z k is denoted by R in Theorem C.1. So R in this case is ||z k || S∞ ≤ R = 4νrL mn + 1 m Further σ is defined by max E(z k z * k ), E(z * k z k )]as per Theorem C.1. We have,EQUATIONNow we would calculate the bound for the above probability, In this scenario, since L = n(n-1), we can bound (σ 2 + Rϵ 3 ) by the following inequality in the second,(σ 2 + Rϵ 3 ) = L m 4νr n + 4νrL mn + 1 m ϵ 3 = 4νrL nm 1 + ϵ 3 + ϵ 3m ≤ 8νrL nm + ϵ 3m ≤ 9νrL nmWe use ϵ 3 is less than 1 in the first inequality and the the second summand is less than νrL nm in the second inequality.Here further we used that ϵ ≤ 1 2 and ν ≥ 1, r ≥ 1 and bound the probability as follows,EQUATIONSo, given that m ≥ 49νnr ϵ 2 log n holds true, we can arrive at the first inequality above, further, since L = n(n-1) 2 so we can say, L ≤ n 2 /2 hence we get the second inequality. We can bound this further by,2n exp -ϵ 2 2 σ 2 + Rϵ 3 ≤ 2n exp -98 log n 18 ≤ 2n exp -4 log n ≤ 2n exp log 1 n 2 = 2 n 3Hence we can conclude from here that || m k=1 z k || S∞ ≤ ϵ) holds with a probability of 1 -2 n 3 . Now we can derive the bound P T0 Q * Ω P T0 -P T0 S∞ ≤ ϵ in the statement of Theorem 4.5 by the following argument,m k=1 z k s∞ = L ℓ=1 ( L m P T0 w α ℓ v ⊤ α ℓ P T0 - 1 m P T0 WV ⊤ P T0 ) s∞ = P T0 R * Ω P T0 -P T0 WV ⊤ P T0 S∞ = P T0 R * Ω P T0 + P T0 w E v ⊤ E P T0 -(P T0 WV ⊤ P T0 + P T0 w E v ⊤ E P T0 ) S∞ = P T0 Q * Ω P T0 -P T0 S∞ ≤ ϵWe get the get implication by adjusting the equation with addition and subtraction of P T0 w E v ⊤ E P T0 because P T0 WV ⊤ P T0 + P T0 w E v ⊤ E P T0 = P T0 W Ṽ⊤ P T0 and W Ṽ⊤ P T0 = P T0 So, we have proved the assertion of Theorem 4.5 is true with a probability of 1 -2 n 3 .Following Proposition 3.1 in [RXH11] , if we change the sampling model to sampling without replacement we should have the same bound with same failure probability.In this section we provide the proof of the local convergence theorem as stated in Theorem 4.3.Lemma D.1. Let 0 < ϵ ≤ 1 2 , let X 0 ∈ S n be a ν-incoherent matrix. Let P T0 : S n → S n be the projection operator associated to T 0 . Then assume that the following three conditions hold:(a) For Q Ω : S n → S n be defined as in (13) from m independent uniformly sampled locations, we have :∥ Q Ω ∥ S∞ ≤ 20L log n m + 1 . (b)The tangent space T 0 = T X 0 onto the rank-r manifold T 0 = T X 0 M r fulfills :EQUATION(c) The spectral norm distance between X and X 0 fulfills:X -X 0 S∞ ≤ ϵ 20L log n m + 1 σ r (X 0 ) (36)Then the tangent space T = T X onto the rank-r manifold at X fulfills:EQUATIONLemma D.1 is a literal adaptation of [KV21, Lemma B.3 ] to the operator Q * Ω , which is why we omit its proof. Similarly, we can use Lemmas B.8 and B.9 of [KV21] in the same way for our proofs.Lemma D.2. Let X 0 ∈ S n be a matrix of rank r that is ν-incoherent, and let Ω = (i ℓ , j ℓ ) m ℓ=1 be a random index set of cardinality |Ω| = m that is sampled uniformly without replacement, or, alternatively, sampled independently with replacement. There exists constants C, C, C 1 such that ifm ≥ Cνrn log n (38)then, with probability at least 1 -2 n 2 , the following holds: For each matrix X (k) ∈ S n fulfillingEQUATION)it follows that the projection P T k : S n → S n onto the tangent space T k := T Tr(X (k) ) M r satisfiesP T k Q * Ω P T k -P T k S∞ ≤ 2 5 ,and furthermore,∥η∥ F ≤ CL log n m ∥P T ⊥ k (η)∥ F . for each matrix η ∈ ker Q Ω in the null space of the operator Q Ω : S n → S n .Proof of Lemma D.2. Assume that there are m locations Ω = (i ℓ , j ℓ ) m ℓ=1 in n × n sampled independently uniformly with replacement, where m fulfills (38) with C := 49/ε 2 and ε = 0.1. By Lemma C.2, it follows that the corresponding operator Q Ω from Lemma C.2 fulfillsEQUATION)on an event called e Ω , which occurs with a probability of at least 1 -2 n , and by Theorem 4.5, the tangent space T 0 = T X 0 M r corresponding to the µ 0 -incoherent rank-r matrix X 0 fulfillsP T0 Q * Ω P T0 -P T0 S∞ ≤ εon an event called e Ω,T0 , which occurs with a probability of at leastEQUATION)it follows by Lemma D.1 that on the event E Ω ∩ E Ω,T0 , the tangent space T k := X (k) onto the rank-r manifold at X (k) fulfillsEQUATION)Next, we claim that on the eventEQUATION)for any for each matrix η ∈ ker Q Ω in the null space of the operatorQ Ω : S n → S n . Let η ∈ ker Q Ω .Then∥P T k (η)∥ 2 F = ⟨P T k (η), P T k (η)⟩ = P T k (η), P T k Q * Ω P T k (η) + P T k (η), P T k (η) -P T k Q * Ω P T k (η) ≤ P T k (η), P T k Q * Ω P T k (η) + ∥P T k (η)∥ F P T k -P T k Q * Ω P T k S∞ ∥P T k (η)∥ F ≤ P k (η), P T k Q * Ω P T k (η) + 4ϵ ∥P T k (η)∥ 2 FUsing (42) in the last inequality, implies thatEQUATIONDividing by ∥P T k (η)∥ F on both sides we get,∥P T k (η)∥ F ≤ 1 1 -4ϵ ∥P T k ∥ S∞ ∥ Q Ω P T k (η)∥ F ≤ 2∥ Q Ω P T k (η)∥ FSince spectral norm of P T k . is bounded by 1. Furthermore, we used that ϵ ≤ 1 10 in the last inequality. Since η ∈ ker Q Ω , it holds that0 = ∥ Q Ω (η)∥ F = Q Ω P T k (η) + P T ⊥ k (η) F ≥ ∥ Q Ω P T k (η)∥ F -∥ Q Ω P T ⊥ k (η)∥ F so that ∥ Q Ω P T k (η)∥ F ≤ ∥ Q Ω P T ⊥ k (η)∥ F ≤ 20L log n m + 1 ∥P T ⊥ k (η)∥ F, where we used (40) in the last inequality. Inserting this above, we obtainEQUATION)for the constant C defined by,C := √ 5 • 21Moreover, we observe that for C 1 := 1 20 C where C 1 is the constant of (39), it holds thatξ ≤ 1 10 • 20L log n m = C 1 √ m L √ log nimplying that the two statements of Lemma D.2 are satisfied on the event E Ω ∩ E Ω,T0 if (39) holds. By the above mentioned probability bounds and a union bound, E Ω ∩ E Ω,T0 occurs with a probability of at least 1 -2n -2 , finishing the proof for the sampling with replacement model. By the argument of Proposition 3 of [Rec11] , the result extends to the model of sampling locations drawn uniformly at random without replacement, with the same probability bound. This concludes the proof of Lemma D.2.The following lemma will also play a role in the proof of Theorem 4.3.Lemma D.3. Let C, C, C 1 be the constants of Lemma D.2 and µ 0 be the incoherence factor of a rank-r matrixEQUATIONProof. First, we compute that∥P T ⊥ k (η (k) )∥ F ≤ ∥P T ⊥ k (X (k) )∥ F + ∥P T ⊥ k (X 0 )∥ F ≤ d i=r+1 σ 2 i (X (k) ) + U (k) ⊥ U (k) * ⊥ X 0 V (k) ⊥ V (k) * ⊥ F ≤ √ n -rσ r+1 (X (k) ) + ∥U (k) * ⊥ U 0 ∥ S∞ ∥Σ 0 ∥ F ∥V * 0 V (k) ⊥ ∥ S∞ ≤ √ n -rσ r+1 (X (k) ) + 2∥η (k) ∥ 2 S∞ (1 -ζ) 2 σ 2 r (X 0 ) √ rσ 1 (X 0 ) = √ n -rσ r+1 (X (k) ) + 2∥η (k) ∥ 2 S∞ (1 -ζ) 2 σ r (X 0 ) √ rκ,where 0 < ζ < 1 such that ∥X (k) -X 0 ∥ S∞ ≤ ζσ r (X 0 ), using Lemma D.4 twice in the fourth inequality and ∥AB∥ F ≤ ∥A∥ S∞ ∥B∥ F all matrices A and B, referring to the notations of lemmaEQUATION⊥ and V⊥ . Using Lemma D.2 for η (k) = X (k) -X 0 , we obtain on the event on which the statement of Lemma D.2 holds that∥η (k) ∥ S∞ ≤ ∥η (k) ∥ F ≤ CL log n m ∥P T ⊥ k (η (k) )∥ F ≤ CL log n m √ n -rσ r+1 (X (k) ) + 8 √ rκ∥η (k) ∥ 2 S∞ σ r (X 0 ) ≤ CL log n m √ n -rσ r+1 (X (k) ) + 8 • 10C 1 √ m √ rκ 8L √ log n √ rκ ∥η (k) ∥ S∞ since C 1 = 1 20 C , after rearranging we get, 1 - 1 2 ∥η (k) ∥ S∞ ≤ CL log n m ( √ n -r)σ r+1 (X (k) )which implies the statement of this lemma.Lemma D.4 (Wedin's bound [Ste06] ). Let X and X be two matrices of the same size and their singular value decompositionsX = (U U ⊥ ) Σ 0 0 Σ ⊥ V * V * ⊥ and X = U U ⊥ Σ 0 0 Σ ⊥ V * V * ⊥ ,where the submatrices have the sizes of corresponding dimensions. Suppose that δ, α satisfying 0 < δ ≤ α are such that α ≤ σ min (Σ) and EQUATIONEQUATIONwithEQUATIONholds with a probability of at least 1 -2n -2 . Then, by lemma B.9 of [KV21] ,EQUATIONWe denote the event that this is fulfilled by E. Furthermore, on this event, if ξ ≤ 1/2 in (48) and denoting the condition number by κ = σ 1 (X 0 )/σ r (X 0 ), it follows from ,lemma B.8 of [KV21] , that∥X (k+1) -X 0 ∥ S∞ ≤ C 2 L log n m 4σ r (X 0 ) -1 ϵ 2 k + 4ϵ k ∥η (k) ∥ S∞ κ + 2∥η (k) ∥ 2 S∞ κ Furthermore, if X (k)r ∈ R n×n denotes the best rank-r approximation of X (k) in any unitarily invariant norm, we estimate thatϵ k ≤ σ r+1 (X (k) ) = ∥X (k) -X (k) r ∥ S∞ ≤ ∥X (k) -X 0 ∥ S∞ = ∥η (k)∥ S∞ , Inserting these two bounds into (50), we obtain∥η (k+1) ∥ S∞ = ∥X (k+1) -X 0 ∥ S∞ = C 2 L log n m 4σ r (X 0 ) -1 (1 + 6κ) ∥η (k) ∥ 2 S∞ .Finally, if, additionally, (48) is satisfied forEQUATIONwe conclude that ∥η (k+1) ∥ S∞ < ∥η (k) ∥ S∞ and also, we observe a quadratic decay in the spectral error such that∥η (k+1) ∥ S∞ ≤ µ∥η (k) ∥ 2 S∞ with a constant µ = m C 2 L log n 1 4(1+6κ)σr(X 0 ). This shows condition b of Theorem 4.3. To show the remaining statement, we can use Lemma D.3 to show that if X (k) is close enough to X 0 , we can ensure that the (r + 1)-st singular value σ r+1 (X (k) ) of the current iterate is strictly decreasing. More precisely, assume now the stricter assumption ofEQUATIONIn fact, if ξ fulfills (49) and (51), we can conclude that on the event E,σ r+1 (X (k+1) ) ≤ ∥η (k+1) ∥ S∞ ≤ ∥η (k+1) ∥ S∞ ≤ C 2 L log n m 4σ r (X 0 ) -1 (1 + 6κ) ∥η (k) ∥ S∞ • ∥η (k) ∥ S∞ < C 2 L log n m 4σ r (X 0 ) -1 (1 + 6κ) √ m √ r 12 CL(log n) √ n -r • ξσ r (X 0 )2 CL log n m ( √ n -r)σ r+1 (X (k) ) ≤ σ r+1 (X (k) )using Lemma D.3 for one factor ∥η (k) ∥ S∞ and (52) for the other factor ∥η (k) ∥ S∞ in the third inequality, and (51) in the last inequality. Taking the update rule (7) for the smoothing parameter into account, this implies that ϵ k+1 = σ r+1 (X (k+1) ), which ensures that the first statement of Theorem 4.3 is fulfilled likewise for iteration k + 1. By induction, this implies that X (k+ℓ) ℓ→∞ ---→ X 0 , which finishes the proof of Theorem 4.3.As a continuation of Section 5, we provide the error analysis for the 1BPM protein data whose corresponding recovery visualizations are found in Figures 7a and 7b . The figure Figure 7a shows the Procrustes distance between the recovered matrix from the samples provided and the ground-truth for both datasets. Similar to the phase transition diagram of the Gaussian data, we observe the probability of success observed over these 24 instances in fig. 7b . Figure 7 Next, we evaluate the performance of MatrixIRLS (Algorithm 1), on the US cities dataset [UU20] in comparison to the aforementioned algorithms. In this setup, we are given m = |Ω| Euclidean distances (λ iλ j ) 2 + (ϕ iϕ j ) 2 between vectors containing longitude and latitude values λ i , λ j and ϕ i , ϕ j of n = 2920 cities in the United States, whose squares can be arranged in an incomplete distance matrix D ∈ S n that serves as an input for the reconstruction algorithms. Like the previous setup of Section 5, the set Ω ⊂ I of point index pairs is sampled uniformly at random, and we consider different choices of m parameterized by the oversampling factor ρ. Here, for the US cities the rank of the input matrix is 2. In Figure 8a we observe a similar pattern to that of Figure 7a , for the US cities data. The success probability Figure 8b also shows a similar pattern to that Figure 4 . We visualize the reconstruction of the US cities data as shown in Figure 9 . 9: Visualizations of the recovery of US cities data by MatrixIRLS for oversampling rates 1.5, 2.5, 3.5 in 9a,9b,9c respectively In this paper, we have adapted the authors' codes of ScaledSGD [ZCZ22] , ALM [TL18] and Ma-trixIRLS [KV21] from their respective githubs. The code of RieEDG [SLCT23] [KV21], has been obtained from the author through personal communication.• MatrixIRLS: The algorithm is configured with the following options:The input parameter for MatrixIRLS includes the number of outer IRLS iterations N 0 , the number of inner iterations N 0 inner , the tolerance, which is the stopping criteria for the algorithm tol 0 inner . For large datasets like the UScities data and Protein data N = 400, although the algorithm converges within 120 iterations. We run the experiments with N 0 inner = 2000 and tol inner = 10 -10 . In a smaller setup of the synthetic data we perform the experiments with n = 500 points with same parameters. although convergence of MatrixIRLS is observed much faster. To put more emphasis on the per iteration error, we study the experiment on per-iterate analysis with more iterations.• Augmented Lagrangian Method: For the Augmented Lagrangian method(ALM) we have same number iterations which is parametrized as N 0 f irstorder for the Barzilai-Borwein gradient method in the code. There are 3 stopping criteria for the BB gradient method which are parametrised by xtol,f tol and gtol on the iterate, functional value of the iterate and the gradient of the function respectively. As mentioned in [TL18] , the relative change in Energy is the stopping criterion for the algorithm and the tolerance for that measure is set at 10 -10 for our experiments. Similar to the above setup, for larger dataset, more iterations are observed so that we can achieve a fair comparison between all the algorithms. • Scaled SGD: This method uses the learning rate and the number of iterations N 0 f irstorder as the parameters. For all the experiments the value of N 0 f irstorder are kept constant across the methods. however, based on the dimension of the input which is n, the iterations are modified to achieve a fair comparison. The learning rate is set at 0.2, which is based on the respective paper [ZCZ22] .• Riemmanian Gardient: This method has parameters, number of iterations which is same N 0 f irstorder , the thresholding tolerance which has same value as the tol inner .The success of proabbility in the experiments of Section 5, is with respect to the Procrustes distance.The Procrustes analysis uses similarity transformations like scaling, rotation and maps into a common reference frame [ADSVF23] . If we consider our problem, this distance is considered with the ground truth as the reference frame and is defined asmin Q∈R r×r ,t∈R r ∥Q • (P 0 + t1 ⊤ ) -P rec ∥ F subject to Q ⊤ Q = I E.. The total degrees of freedom in D is r. Hence, the total degrees of freedom of A is: nr + r -(r + r(r + 1) 2 ) = nr -r(r -1) 2 .For a symmetric matrix X ∈ S n , we write its eigendecomposition (with in magnitude decaying eigenvalues) such thatEQUATIONwhere U ∈ R n×r , U ⊥ ∈ R n×(n-r) are matrices with orthonormal columns and Λ = diag(γσ 1 , . . . γσ r ) and Λ ⊥ := dg(γσ r+1 , . . . γσ d ) diagonal matrices with entries λ i = σ i γ i where σ i is the i-th singular value of X and γ i ∈ {±1} a sign. Furthermore, we denote by T r (X) the best rank-r approximation of X, which can be written such that T r (X) := arg minEQUATIONwhere ∥ • ∥ can be any unitarily invariant norm, due to the Eckardt-Young-Mirsky theorem [Mir60], using the established notation.The computational complexity of Algorithm 1 is dominated by the solution of the weighted least squares problem (5) and the computational of spectral information used in the update of the smoothing parameter ϵ k of (7) and the weight operator update (7) (see also Definition 3.1).For solving (5), we consider detailed in Appendix F.1, based on solving a linear system with the dimensionality of tangent space T Tr k (X (k) ) M r k . The tangent space formulation can be derived from (6) via the Sherman-Morrison-Woodbury [Woo50] formula using the weight operator structure.In order to update the smoothing parameter and the weight operator we compute the SVD of the iterate by approximating top singular vectors and corresponding values following Randomized Block Krylov Methods [MM15] . This computation takes O(((m + nr)r log n).Let M r := {X ∈ S n : rank(X) = r} the manifold of symmetric rank-r matrices and let X ∈ S n be as in (53). In this case, given r ∈ N, we can write the tangent space of M r at T r (X) asEQUATIONsee also [Van13] , [Bou20, Chapter 7.5].If P T : S n → S n is the orthogonal projection operator that projects symmetric matrices onto T as used in Theorem 4.5, we note that P T (X) = UU ⊤ XUU ⊤ + UU ⊤ X(I -UU ⊤ ) + (I -UU ⊤ )XUU ⊤ , which can be decomposed such that P T (X) = P T P * T (X), where the action of P T : R r(n+r) → S n can be described asP T (γ) := UΓ 1 U ⊤ + UΓ ⊤ 2 + Γ 2 U ⊤ with Γ 1 ∈ S rbeing the result of an (r 2 × 1) to (r × r) reshaping of the first r 2 coordinates of γ and Γ 2 ∈ R n×r the result of an (rn × 1) to (n × r) reshaping of the remaining coordinates of γ. Further, we used the notation of the adjoint operator P * T : S n → R r(n+r) of P T which maps X ∈ S n to the vectorization γ of [Γ 1 , Γ 2 ] := U ⊤ XU, (I -UU ⊤ )XU .Due to the fact that the weight operator W k = W X (k) ,ϵ k of Definition 3.1 associated to the iterate X (k) ∈ S n and the smoothing parameter ϵ k > 0 is self-adjoint and invertible, we can write its inverse asEQUATIONwhereD -1 k ∈ R r k (n+r k )×r k (n+r k) is a diagonal matrix with entries max(σ i (X (k) ), ϵ k ) max(σ j (X (k) ), ϵ k ), where either i or j are smaller or equal than r k , see also [KV21, Eq. ( 12)], and T k = T Tr k (X (k) ) M r k is the tangent space at T r k (X (k) ) onto the rank-r k manifold.Recall from (7) that the the defining equation for the next iterate X (k+1) given W k isEQUATIONUsing (56), we see thatEQUATIONand, using the Sherman-Morrison-Woodbury [Woo50] formula, we have thatEQUATIONwith linear system matrix M :=ϵ 2 k C -1 + P * T k A * (AA * ) -1 AP T k , noting that C := D -1 k -ϵ 2 k I is invertible since (D -1 k ) ii > ϵ 2k for all diagonal indices. The observation of (59) can be turned an efficient implementation of the weighted least squares computing X (k+1) , which is presented in Algorithm 2. It can be shown that Algorithm 2 indeed computes X (k+1) implicitly, cf. Lemma F.1. We omit its proof, which follows the proof of [KV21, Lemma A.1].Algorithm 2 Tangent space implementation of weighted least squares step in Algorithm 1Input: Set Ω, observation distances D Ω = (d ij ) (i,j)∈Ω , singular vectors U ∈ R d1×r k , V (k) ∈ R n×r k and singular values σ (k) 1 , . . . , σ (k) r k , smoothing parameter ϵ k , projection γ (0) k = P * T k P T k-1 (γ k-1 ) ∈ R r k (n+r k ) of solution γ k-1 ∈ R r k-1 (n+r k-1) of linear system (60) for previous iteration k -1.EQUATION) for ∆γ k ∈ S k by the conjugate gradient method [HS52, Meu06] , whereM = ϵ 2 k D -1 k -ϵ 2 k I -1 + P * T k A * (AA * ) -1 AP T k . 3: Compute γ k = γ (0) k + ∆γ k . 4: Compute residual r k+1 := (AA * ) -1 (y -AP T k (γ k )) ∈ R m . Output: r k+1 ∈ R m and γ k ∈ R r k (n+r k ) . Lemma F.1. If r k+1 ∈ R m and γ k ∈ R r k (n+r k )is the output of Algorithm 2, then X (k+1) as in step (6) of Algorithm 1 fulfillsX (k+1) = A * (r k+1 ) + P T k (γ k ) = A * (r k+1 ) + UΓ 1 U ⊤ + UΓ ⊤ 2 + Γ 2 U ⊤ , where Γ 1 ∈ R r k ×r k is the matricization of the first r 2k elements of γ k and Γ 2 ∈ R n×r k the matricization of the reamining elements.Algorithm 3 Implementation of U ⊤ A * : R m+n → R r×n Input: Input vector y = [y 1 , y 2 , ....., y m+n ], Ω = (i ℓ , j ℓ ) m ℓ=1 ⊂ I be a multiset of double indices. S 1 = m ℓ=1:(i ℓ ,j)∈Ω for some j y ℓ e i ℓ e T i ℓ S 2 = m ℓ=1:(i,j ℓ )∈Ω for some i y ℓ e j ℓ e T j ℓ S 3 = m ℓ=1:(i ℓ ,j ℓ )∈Ω y ℓ e i ℓ e T j ℓ S = S 1 + S 2 -S 3 -S T 3 . saved as a sparse matrix.V 1 = U ⊤ S ▷ 4rm flops. V 2 = 1 2 U ⊤ • [y m+1 , ..., y m+n ] ⊤ • 1 ⊤ ▷ rn flops. V 3 = 1 2 U ⊤ • 1 • [y m+1 , ..., y m+n ] ▷ 2rn flops V = V 1 + V 2 + V 3 Output: V Total of rm + 3rn flops.Algorithm 4 Implementation of AA * : R m+n → R m+n Input: Ω = (i ℓ , j ℓ ) m ℓ=1 ⊂ I be a multiset of double indices fulfilling m < L = n(n -1)/2 that are sampled independently with replacement, the vector y = [y 1 , y 2 , ....., y m+n ].Load S from Algorithm 3 which computes A * (y).T 1 = (e ⊤ i ℓ Se i ℓ ) m ℓ=1:(i ℓ ,j)∈Ω for some j T 2 = (e ⊤ j ℓ Se j ℓ ) m ℓ=1:(i,j ℓ )∈Ω for some i . EQUATION∈ R n×r . A 1 = V ∈ R r×n from Algorithm 3. ▷ O(r(m + n)) flops Γ 1 = A 1 U ∈ R r×r . ▷ O(r 2 n) flops M = Γ 1 U ⊤ ∈ R r×n . ▷ O(r 2 n) flops Γ 2 = (A ⊤ 1 -M ⊤ ) Output: {Γ 1 , Γ 2 }.Lemma F.1 shows that an iterate X (k+1) can be represented via only m + r k (n + r k ) parameters. In the remainder of this discussion, we will assume that r k = r, which is the case in most cases in practice if the rank estimate r of Algorithm 1 is chosen as r = r.To quantify the computational cost of Algorithm 2, we assume a fixed number N 0 inner of CG iterations solving (60). When applying the system matrix M via matrix-vector multplication, we observe that its first summand ϵ 2 k D -1 kϵ 2 k I -1 is diagonal and thus results in O(r(n + r)) flops per CG iteration.To quantify the matrix-vector multiplication cost of its second summand P * T k A * (AA * ) -1 AP T k , we define below algorithms that efficiently implement the application of the operators P * T A * : R m+n → R r(n+r) (Algorithm 5), AA * : R m+n → R m+n (Algorithm 4) and AP T : R r(n+r) → R m+n (Algorithm 6). As an auxiliary function, we also provide an implementation of U ⊤ A * : R m+n → R r×n (Algorithm 3) below. The application of (AA * ) -1 : R m+n → R m+n can be achieved by an inexact iterative solver that applies Algorithm 4 a fixed number Ñ of iterations, which costs O((m + n) Ñ ). Since the time complexity of Algorithm 5 and Algorithm 6 are O(r 2 n + rm), we obtain one matrix vector multiplication with M from Algorithm 2in a time complexity of O(r 2 n + rm + (m + n) Ñ ).Since N 0 inner CG iterations are used, we obtain a total time complexity of O(N 0 inner (r 2 n+rm+(m+ n) Ñ )) for Algorithm 2. In our experiments, we observe that a small number Ñ = 10 of iterations for solving the system associated to (AA * ) -1 is sufficient to obtain high-accuracy solutions. This breakdown of the computational costs of each algorithm is shown in Algorithm 3 to Algorithm 6.Algorithm 6 Implementation of AP T : R r(n+r) → R m+nInput: {Γ 1 , Γ 2 } , index set Ω, left singular vectors U ∈ R n×r . M = Γ 1 U ⊤ ∈ R r×n . ▷ O(r 2 n) flops from Algorithm 5 z = ( r k=1 U (j ℓ ,k) M (k,j ℓ ) ) (i,j ℓ )∈Ω for some i ▷ O(mr) flops z = z + ( r k=1 U (i ℓ ,k) M (k,i ℓ ) ) (i ℓ ,j)∈Ω for some j ▷ O(mr) flops z = z -( r k=1 U (j ℓ ,k) M (k,i ℓ ) ) (i ℓ ,j ℓ )∈Ω ▷ O(mr) flops z = z -( r k=1 U (i ℓ ,k) (M (k,j ℓ )) (i ℓ ,j ℓ )∈Ω ▷ O(mr) flops α = ( r k=1 U (i ℓ ,k) (Γ ⊤ 2 ) (k,i ℓ ) ) (i ℓ ,j)∈Ω for some j ▷ O(mr) flops α = α + ( r k=1 U (j ℓ ,k) (Γ ⊤ 2 ) (k,j ℓ ) ) (i,j ℓ )∈Ω for some i ▷ O(mr) flops α = α -( r k=1 U (j ℓ ,k) (Γ ⊤ 2 ) (k,i ℓ ) ) (i ℓ ,j ℓ )∈Ω ▷ O(mr) flops α = α -( r k=1 (U) (i ℓ ,k) (Γ ⊤ 2 ) k,j ℓ ) (i ℓ ,j ℓ )∈Ω ▷ O(mr) flops ζ 1 = z + 2α . c 1 = ( n i=1 U (i,k) ) r k=1 ▷ O(nr) flops γ 1 = c 1 M ▷ O(mr) flops γ 2 = c 1 Γ ⊤ 2 ▷ O(mr) flops c 2 = ( n i=1 (Γ 2 ) (i ℓ ,k) ) r k=1 ▷ O(nr) flops γ 3 = c 2 U ⊤ ▷ O(mr) flops ζ 2 = γ 1 + γ 2 + γ 3 Output: [ζ 1 ; ζ 2 ]• The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions). • The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.) • The assumptions made should be given (e.g., Normally distributed errors). • It should be clear whether the error bar is the standard deviation or the standard error of the mean. Guidelines:• The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.• If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction).For example, for a rank r and oversampling factor ρ, if out of theinstances, in 12 such cases, the Procrustes distance based error is less than the threshold, then the probability of success is 1 2 . 


Decision: The submission studies the recovery of a low-rank Euclidean distance matrix from a subset of its entries. This problem has applications in sensor positioning, protein folding, and related areas. The paper studies an iterative reweighed least squares algorithm, which solves a sequence of weighted least squares completion problems. It proves that this method converges locally, with a linear rate. In particular, if the observed distances are sampled uniformly at random, and the number of observations is at least n r log n times the coherence of the target matrix, then IRLS converges linearly within a certain schatten-infinity neighborhood of the truth. In experiments, the IRLS approach shows advantages in handling ill-conditioned data — a challenging case for nonconvex factorization-type methods. 

Reviewers produced a generally positive evaluation of the paper: while there has been quite a bit of theoretical work on nonconvex low-rank recovery, including analyses of both factorization and iterative least squares methods, the paper has novelty within the literature on nonconvex *euclidean distance matrix* reconstruction. The paper is crisply written. The main limitation of the theory is that it does not comprise a complete recovery result, since convergence is local in nature. One device for achieving a recovery result is to couple the local analysis with spectral initialization; doing so would help to contextualize the convergence guarantee (i.e., is the radius large enough to ensure near optimal sample complexity for this hybrid method). Overall, the problem makes a solid contribution to the literature on nonconvex low-rank recovery.

Reviews: review: paper_summary: paper_summary: The paper examines the problem of recovering point locations based on incomplete pairwise Euclidean distance information. The problem is cast as a low-rank matrix recovery problem, and an algorithm based on iteratively-reweighted least squares is proposed. The authors derive a lower bound for the minimum number of required samples for reconstruction, and show that their algorithm has quadratic convergence. Experimental results show that their method has comparable results to other state-of-the-art methods for most data, but is also uniquely robust to ill-conditioned cases where other methods fail to converge.
strength_weakness: strength_weakness: * Example data is somewhat small-scale and synthetic, but appears sufficient for comparison with other recovery methods. 
* While this is not a weakness of the paper, my own theoretical background is insufficient to verify some of the proofs, such as the construction of the dual basis followed by the establishment of the restricted isometry property condition.
questions: questions: * While number of steps/time to convergence is shown for ill-conditioned matrix experiments, they are not shown for real data experiments. How much faster does the proposed MatrixIRLS derivative converge compared to ALM? Are they comparable?
* Are there applicable downstream machine learning/deep learning applications that require this kind of point recovery, where improved results could be further demonstrated?
limitations: limitations: The authors address the limitations of their work in detail in the paper, and there are no societal impacts to be addressed.
paper_summary: paper_summary: This paper tackles the problem of reconstructing geometric configurations from minimal Euclidean distance samples. The authors propose an innovative algorithm based on iteratively reweighted least squares (IRLS) within a non-convex optimization framework. They provide a local convergence guarantee for the method, ensuring accurate geometry recovery from a minimal number of random distance samples. Additionally, the paper establishes a restricted isometry property (RIP) relevant to other non-convex approaches. Numerical experiments demonstrate the algorithm’s efficiency and robustness compared to state-of-the-art methods.
strength_weakness: strength_weakness: 1. The paper compares the proposed method with a few state-of-the-art algorithms. Including more diverse comparisons, particularly with newer methods covered in the existing literature, would provide a more comprehensive evaluation.

2. The experiments could be enhanced by adding more parameter analysis or robustness analysis. For example, testing the performance under different parameters (n, r, k) on the ill-conditioned data will further validate the effectiveness of the proposed method.

3. The paper provides a local convergence guarantee, but a more detailed discussion on the global convergence properties of the algorithm would be beneficial.
questions: questions: 1.How sensitive is the proposed algorithm to the choice of parameters, such as rank estimate in the IRLS framework? Providing guidelines or heuristics for parameter selection would be helpful.

2.In the additional experiments presented in Appendix E.1, why does the MatrixIRLS method exhibit a large standard deviation when rho is in the range of 2.5 to 3.5? This suggests that the method is not very stable or robust across all settings. Considering the success probability, ALM performs better than MatrixIRLS. Could you explain the potential reasons for this discrepancy?

3.The time complexity of the MatrixIRLS algorithm appears to be missing a component in Step 7, which involves singular value decomposition with a complexity of O(n^3). Given this, what is the scalability of the MatrixIRLS method, such as the maximum problem size that can be solved within one hour?

4.Will the results of the MatrixIRLS method converge to the same solution as MDS when there is no incomplete information?

*Small suggestions:*

1.The linewidth and point size is too large in Figs. 2, 6(b), and 7(b).

2.No captions in Figs. 6 and 7.

3.Mixed use of fig. and Figure.

4.Some typos: 

Line 52: should satisfy also satsify

Line 813: figure Figure
limitations: limitations: It would be beneficial to include the scalability limitations, such as the maximum problem size that can be solved.
paper_summary: paper_summary: The paper studies the under-constrained problem to reconstruct n ordered points in Euclidean space R^r, which are given by fewer than all n(n-1)/2 pairwise distances.
strength_weakness: strength_weakness: Different from the initial Problem 1, the actual considered problem is the rank minimization in (1) after line 54. Does the paper explain anywhere why it makes sense to minimize this integer rank of a Gram matrix even if we forget about initial points. Another weakness is the practical impossibility of solving this rank minimization, which is claimed in line 54 to be NP-hard. It would be useful to re-phrase a problem in another way that allows a polynomial solution. 

More importantly, the rank minimization in (1) and algorithm MatrixILRS appeared under the same name but with slightly different notations in the earlier paper [KV21]. So the problem and algorithm in the submitted paper look incremental updates of the past work, which even uses the same names, includes similar theorems, example structures and pseudo-codes.
questions: questions: For a simple example, if n=4 points in the plane are given by 4 distances that are edge lengths of a quadrilateral on these 4 points, then what's the meaning of this rank minimization in this case? Can Problem 1 be explicitly solved in some partial cases, e.g. for 5 distances on 4 points in the plane? 

In Def 3.2, is w_alpha a matrix n x n? In lines 177-178, how is measurement operator <w_alpha_l,X> computed? Are angle brackets missed for l>m in line 178?  

The concept of a nu-incoherent matrix appears in Theorems 4.2, 4.5 without earlier explanations. What are these matrices?

The words "coherence factor" appear for the first time in line 80 without defining this factor nu. Also, nu appears earlier in line 66, still without explanations. Is it the same "coherence" appearing much later in Definition 4.1?
 
The word "probability" appears for the first time in Theorem 4.2, later again in Theorem 4.5. How is the probability space defined to make the convergence with "high probability" meaningful?

In Theorem 4.3, can mu be large enough so that the resulting inequality becomes equality and hence no convergence is guaranteed? What are the values of mu in the experiments form section 5?

Lines 251-252 say that "an adaptation of the proof of [KV21] to the EDG setting is not possible due to non-orthogonality of the basis {w_alpha}_alpha of Definition 3.2". What specific step in the past proof requires this orthogonality?

Lines 857-858 say that "The Procrustes analysis uses similarity transformations like scaling, rotation and maps into a common reference frame [ADSVF23]". However the minimization after line 859 includes only translations and orthogonal maps. Was uniform scaling missed?  Why is the minimum achieved over infinitely many possible translations? Does this Procrustes distance satisfy the metric axioms and what is its computational complexity in the number of points?

The biggest concern is line 82 saying that "from which the geometry of the points P is trivially recovered". Main Algorithm 1 (MatrixIRLS) outputs a matrix X^(k) through some iterative optimization. Why is this final n x n matrix realizable by n points in a Euclidean space R^r? For instance, we cannot arbitrarily choose 6 pairwise distances between 4 points in the plane. Even if these distances are positive and satisfy triangle inequalities, they should also satisfy a complicated polynomial equation saying that the tetrahedron on these 4 points has volume 0. For n>4 points, there are many more such equations, so a random symmetric matrix of positive numbers is realized by points with probability 0.
limitations: limitations: This submission seems an incremental improvement of the past work [KV21] with little comparison between the two papers.
paper_summary: paper_summary: This paper proposed a method called MatrixIRLS to solve the EDG problem, an important problem in various fields, especially when pairwise distance is partially observed. The proposed method improves the number of random sampled distances to guarantee local convergence to the ground truth.
strength_weakness: strength_weakness: 1.	There are quite a few typos or undefined notation that make some of the paragraphs hard to read. For example:

  a.	Line 66, $\nu$ is not defined. Although it’s briefly mentioned in line 80, its rigorous definition is on Page 5. I suggest the following change to line 66: $O()$ (where $\nu$ is the coherence factor, see \cref{def:4.1} for details).

  b.	Line 158, $\sigma_i(\cdot)$ is not defined, tho I can guess the meaning.

  c.	Line 165: $X^{(k)}$ should be $X$. 

  d.	Eq (5) gives the update from k to k+1, while Eq (6) gives the rule from k-1 to k.

  e.	Eq(9), at the end, there is a missing comma.

2.	Figure 2 looks great, but it will be even better if we can observe the claimed rate $O(\nu rn\log n)$ empirically.
questions: questions: 1.	There are some vertical space issues, like Line 286 for Section 5. 
2.	Can the authors briefly discuss the main challenges for global convergence?
limitations: limitations: Discussed in the supplementary material. I suspect the file in the supplementary zip file should be the main submitted file.

Comments: title: Thanks
comment: Thank you for your detailed clarification. As my concerns have been addressed, I will maintain my positive score. By the way, I would like to note that the time complexity you mentioned is not linear in $n$, but rather $O(n \log n)$.
I thank the authors for their work on the rebuttal. Based on their answers to my review, which are further positive points for the paper, as well as reading and considering the other reviews and their rebuttals, I see no reason to change my positive rating of the paper.
If you feel that we have addressed your questions and concerns adequately, we would appreciate if you can consider updating your score. We warmly welcome any further questions that you may have. Thank you!
> To summarize: symmetry holds, but positive-definiteness of the result is true with high probability only if there exists a ground truth matrix which is reconstructed, but there is no proof that the optimization result is indeed positive-definite in all reconstruction scenarios.

That is correct. As previously mentioned, in such cases, one can obtain a realizable point configuration by applying a truncated eigendecomposition. This point configuration will not satisfy the distance constraints exactly anymore, but approximately. Comparing such an outcome with the results of reference algorithms such as ALM, RieEDG and ScaledSGD, we would like to point out that the iterates of these algorithms do not satisfy the distance constraints in the first place, unless their objectives converge to 0, which they do not for many challenging problem instances. On the other hand, the satisfaction of the distance constraints by all iterates is built into MatrixIRLS. Thus, applying a truncated eigendecomposition to the outcome of MatrixIRLS is not less justified than using any of the other algorithms.
Thank you for your clarifications. 
To summarize: symmetry holds, but positive-definiteness of the result is true with high probability only if there exists a ground truth matrix which is reconstructed, but there is no proof that the optimization result is indeed positive-definite in all reconstruction scenarios. 
This coincides with my initial understanding.
Thanks for the response.

Your explanation to the rate makes sense to me. It will be very helpful if the authors can discuss this in section 5, to call back the claimed rate in Section 1.1 

Similarly, discussing the challenges for global convergence will make the story more complete.

I raised the score to 6, mainly because the rate is explained clearly to me, which is part of the main contribution claimed in Section 1.1.
Regarding symmetry, we want to note that all iterates $\mathbf{X}^{(k)}$ of Algorithm 1 are symmetric by construction: By (6), the solution of the equality-constrained weighted least squares problem (5) (with $W_{k-1}$ instead of $W_{k}$) can be written as
$$
\mathbf{X}^{(k)} = W_{k-1}^{-1}\mathcal{A}^*(z_{k-1}),
$$
where $z_{k-1}:= (\mathcal{A}W_{k-1}^{-1} \mathcal{A}^*)^{-1}(\mathbf{y})$ is the solution of the linear system $\mathcal{A}W_{k-1}^{-1} \mathcal{A}^* \mathbf{z} = \mathbf{y}$. From the definition of the measurement operator $\mathcal{A}$ (see line 177-178), it becomes clear that the adjoint operator $\mathcal{A}^*$ of $\mathcal{A}$ maps a vector $z \in \mathbb{R}^{m + n}$ to $\mathcal{A}^*(z) = \sum_{\alpha_{\ell} \in \Omega} z_\ell w_{\alpha_\ell} + \sum_{i=1}^n z_{m+i} w_{(i-m,i-m)}$, where the $w_{\alpha}$ are as in (4). Since all $w_{\alpha}$ are symmetric, it holds therefore that also $\mathcal{A}^*(z_{k-1})$ is symmetric as a sum of symmetric matrices. Based on the definition $W_{k-1}: = W_{\mathbf{X}^{(k-1)},\epsilon_{k-1}}$ of $W_{k-1}$ and Definition 3.1, it can be seen that $W_{k-1}^{-1}$ maps a matrix $\mathbf{Z}$ to $\mathbf{U}\Sigma_{\epsilon_{k-1}}\mathbf{U}^*\mathbf{Z}\mathbf{U}\Sigma_{\epsilon_{k-1}}\mathbf{U}^*$, where $\mathbf{U}$ is a left singular vector matrix of $\mathbf{X}^{({k-1})}$, and $\Sigma_{\epsilon_{k-1}}  = \operatorname{diag}(\max(\epsilon_{k-1},\sigma_i(\mathbf{X}^{(k-1)}))$, which is itself symmetric if $\mathbf{Z}$ is symmetric. Thus, each iterate $\mathbf{X}^{(k)}$ of Algorithm 1 is symmetric.

Regarding positive semidefiniteness, it is true in general, iterates $\mathbf{X}^{(k)}$ can have negative eigenvalues. However, we prove in Theorem 4.3 that if Algorithm 1 is initialized close enough to a rank-$r$ ground truth matrix $\mathbf{X}^0$ and if the sample complexity satisfies the bound of line (229) under the mentioned random sampling model, with high probability, $\mathbf{X}^{(k)}$ converges to $\mathbf{X}^0$ (with quadratic rate). For the problems considered, the ground truth $\mathbf{X}^0 \in S_n$ is a realizable rank-$r$ matrix and thus p.s.d., so the iterates $\mathbf{X}^{(k)}$ converge to a psd matrix under these conditions. The proof of this theorem can be found in lines 786 to 809 (but uses also other intermediate lemmas, cf. the supplementary material).
Thank you for the detailed reply.

>As the reviewer points out, not every random symmetric matrix corresponds to a set of points; positive definiteness is required to realize such points. Under the assumptions of our main theorem, the proposed algorithm converges to the true Gram matrix of rank $r$, which is both symmetric and positive definite. We then obtain the configuration of the points through an r-truncated eigendecomposition.

Is it possible to specify the line numbers where the final matrix is proved to be symmetric and positive definite?