Optimization Can Learn Johnson Lindenstrauss Embeddings

Authors: nikos  tsikouras
constantine  caramanis
christos  tzamos
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-15 17:22:35

Abstract: Embeddings play a pivotal role across various disciplines, offering compact representations of complex data structures. Randomized methods like Johnson-Lindenstrauss (JL) provide state-of-the-art and essentially unimprovable theoretical guarantees for achieving such representations. These guarantees are worst-case and in particular, neither the analysis, ${\textit{nor the algorithm}}$, takes into account any potential structural information of the data. The natural question is: must we randomize? Could we instead use an optimization-based approach, working directly with the data? A first answer is no: as we show, the distance-preserving objective of JL has a non-convex landscape over the space of projection matrices, with many bad stationary points. But this is not the final answer. 

We present a novel method motivated by diffusion models, that circumvents this fundamental challenge: rather than performing optimization directly over the space of projection matrices, we use optimization over the larger space of $\textit{random solution samplers}$, gradually reducing the variance of the sampler. We show that by moving through this larger space, our objective converges to a deterministic (zero variance) solution, avoiding bad stationary points. 

This method can also be seen as an optimization-based derandomization approach, and is an idea and method that we believe can be applied to many other problems.

Main paper: 
 introduction: 
 Embeddings are foundational across diverse disciplines, offering compact representations of complex data structures. Algorithms across different domains leverage embeddings to capture nuanced relationships between data points, improving efficiency and effectiveness in processing. Within this field, there are distinct lines of work.Embeddings have been used as dimensionality-reducing tools, while preserving important structures in data. They have been a major research focus for years, and have been a key ingredient in several algorithmic applications such as graph sparsification [Spielman and Srivastava, 2008] , nearestneighbor search [Indyk and Motwani, 1998] , hashing [Dietzfelbinger et al., 1997] and digital research [Schmidt, 2018] . A celebrated result in this area has been the Johnson-Lindenstrauss (JL) lemma Johnson [1984] which shows that a random linear mapping can reduce the dimension of the dataset while approximately preserving the L 2 norm of all points, with high probability. JL embeddings give bounds on the maximum distortion over all the points. Many variants have been developed and studied, as we discuss below in Section 2. For this paper, the salient point is that the algorithms that enjoy theoretical guarantees are for random projections drawn from a distribution over projection matrices. To the best of our knowledge, even among derandomization of JL-type techniques (e.g., Clarkson and Woodruff [2009] , and see Section 2 for further discussion), there are no guarantees for optimization-based algorithms that directly attempt to minimize a distortion objective.On the other hand, in the last decade, embeddings have risen to prominence across various tasks in deep learning. In encoder-decoder architectures, embeddings act as intermediary representations to solve a broad range of challenges in natural language processing and speech processing, [Bengio and Heigold, 2014 , Cho et al., 2014a ,b, Qi et al., 2018 , Rush et al., 2015 , Xiong et al., 2016] . In the domain of face recognition, embeddings encode facial features into concise vectors, enabling precise identification and matching [Liu et al., 2015 , 2017 , Schroff et al., 2015] . Furthermore, contrastive learning techniques utilize embeddings to emphasize disparities between similar and dissimilar instances, thereby augmenting the discriminative capabilities of models [Gao et al., 2021 , Khosla et al., 2020 , Radford et al., 2021] .In contrast to the JL-type results, in the above applications these embeddings are learned using optimization as part of the (pre)training process [Caron et al., 2021 , Oquab et al., 2023 , Press and Wolf, 2016 , Vaswani et al., 2017] . Though their empirical success is remarkable, the non-convex landscape of the optimization process makes obtaining theoretical guarantees a key challenge.Another important related area is the extensive body of literature dedicated to employing optimization techniques for Principal Component Analysis (PCA). PCA aims to find linear embeddings that capture the optimal average distortion in the data, thereby reducing its dimensionality while preserving as much information as possible. The last decade saw significant success of direct matrix optimization in various PCA-like settings [Garber and Hazan, 2015 , Shamir, 2016 , Xu and Li, 2021 , Yi et al., 2016] . Though clearly there are similarities to the JL objective, the difference between guarantees on the maximum perturbation (JL) vs the average perturbation (PCA) are significant, and one reason why the techniques pioneered for PCA have yet to be applied successfully to JL.Nevertheless, the empirical success of optimization in neural networks, and its theoretical success for PCA-type objectives, motivate us to revisit the JL low-distortion embedding task (see Section 3 for the exact definition of the JL guarantee), and ask the question that to the best of our knowledge, has yet to find an answer:Can optimization-based approaches be used to obtain a matrix that satisfies the Johnson-Lindenstrauss guarantee?Answering this question is the main goal of this paper.Our Contributions and a Conceptual Road Map.The main contribution of this paper is in developing a framework that allows a direct (and deterministic) optimization approach to obtain the same JL guarantees as random projection. After all, it is well-known that the JL guarantees given by random projection are not improvable [Larsen and Nelson, 2014 , 2017 , Alon and Klartag, 2017] . Moreover, established derandomization techniques based on conditional expectation and other methods have long been available [Raghavan, 1988 , Engebretsen et al., 2002 , Bhargava and Kosaraju, 2005] .The key conceptual steps on the way to our main result are as follows:Step 1: We first show that the optimization landscape in the ambient space of projection matrices is not favorable, and in particular, any attempt to directly minimize distortion over this space using first or second-order methods, is destined to fail.Theorem (Informal version of Theorem 1). The maximum distortion objective considered as a function in the space of matrices has many suboptimal local minima.Step 2: Given the above negative result, a different approach is required. We draw inspiration from diffusion models and solution samplers [Bello et al., 2016 , Ho et al., 2020] . Rather than optimize in the space of matrices, we optimize in the larger space of (mean, variance) parameters of Gaussian distributions over embedding matrices. Thanks to the original JL theorem, we know an initial choice of parameters that define a solution sampler whose expected distortion is small: zero mean and unit variance. Akin to diffusion, we then seek to sequentially decrease the magnitude of the variance, without increasing the expected distortion of the sampler. Note that the space of matrices is properly contained in this space of samplers, as we can identify a specific (deterministic) projection matrix with a Gaussian distribution with that mean, and zero variance. The challenge is to find a path from our initial JL sampler, to a deterministic sampler (a projection matrix) whose maximum distortion is approximately as good as the expected guarantee of the original JL sampler.We turn this into an optimization problem by creating an objective function in the space of samplers.Our first result demonstrates that if we find a second-order stationary point for this objective function, then we have solved our original problem. Specifically:Theorem (Informal Version of Theorem 2). For data x 1 , . . . , x n ∈ R d and target dimension k all second-order stationary points reachable from the origin for the objective function defined in Equation 4 have zero variance and hence correspond to fixed matrices; moreover, these matrices satisfy the JL guarantee.Step 3: The final step requires proving the tractability of finding a second-order stationary point, i.e., of solving this optimization problem. We do so using a generic deterministic second-order optimization algorithm (see Alg. 1). Thus our result shows that our second-order optimization algorithm successively performs reverse-diffusion-like steps, decreasing the variance without deteriorating the quality of the solution sampler, until it has finally arrived at a deterministic solution.Theorem (Informal Version of Theorem 3). For data x 1 , . . . , x n ∈ R d and target dimension k running Algorithm 1 using Equation 4for poly(n, k, d) steps returns a matrix that satisfies the JL guarantee.Step 4: Finally, we show through simulations that the qualitative and quantitative results of our theory are borne out. 
related work and alternative approaches: 
 The JL lemma is a well-studied result studied in the literature, with several simplifications and extensions of the original proof [Dasgupta and Gupta, 2003 , Frankl and Maehara, 1988 , Kane and Nelson, 2010 , Matoušek, 2008] . It has also been shown that the JL lemma is optimal in terms of the target dimension. The authors in [Larsen and Nelson, 2017, Alon and Klartag, 2017] provide a tight lower bound on the target dimension required by the JL lemma, given a specific distortion, for any random linear mapping. In addition to these theoretical insights, there have been various approaches aimed at efficiently constructing random matrices that satisfy the JL guarantee with high probability.One approach samples each matrix entry independently from a Gaussian distribution [Indyk and Motwani, 1998 ], while others utilize Rademacher and sparse Rademacher distributions [Arriaga and Vempala, 2006, Achlioptas, 2001] . Moreover, generalized sampling methods have shown that any distribution with zero mean, unit variance, and a subgaussian tail can be used [Matoušek, 2008] . The Fast Johnson-Lindenstrauss Transform employs sparse matrices and randomised Walsh-Hadamard transforms for efficiency [Ailon and Chazelle, 2009] . Additionally, the subsampled randomised Hadamard transform, achieves efficient embeddings by combining subsampling with randomised Hadamard transforms, maintaining a high probability of preserving the distances between points [Ailon and Liberty, 2013] .Derandomizing Johnson-Lindenstrauss. Derandomization is a technique for developing deterministic algorithms or algorithms that require fewer random bits, and it has proven to be a powerful theoretical tool [Kabanets, 2002] . There have been numerous efforts to derandomize the JL lemma, with a significant focus on using pseudorandom generators (PRGs) capable of fooling statistical tests [Nisan, 1990] . These constructions aim to achieve reduced seed lengths while satisfying the JL lemma's norm-preserving properties with ±ε distortion and probability of failure δ. For example, the ℓ 2 -streaming algorithm achieves a JL family with seed length O(log d) and with k = O(1/(ε 2 δ)) [Alon et al., 1996] . The authors in [Clarkson and Woodruff, 2009] leveraged the use of scaled random Bernoulli matrices with Ω(log(1/δ))-wise independent entries, resulting in a seed length of O(log(1/δ) log d). Additionally, PRGs that δ-fool degree-2 polynomial threshold functions generate JL families with seed lengths of poly(1/δ) log d [Meka and Zuckerman, 2010] .Other approaches have utilized conditional probabilities and pessimistic estimators, introduced in [Raghavan, 1988] , to derive deterministic algorithms for JL embeddings [Engebretsen et al., 2002, Bhargava and Kosaraju, 2005] .These works differ from ours in several ways. The Nisan pseudorandom generator uses a few random bits. Additionally, the authors in [Bhargava and Kosaraju, 2005] fully derandomize the Rademacher construction by [Achlioptas, 2001] using a combinatorial algorithm that greedily selects the best matrix entries. Even in such a coordinate-wise fashion, using a gradient-descent (continuousoptimization) approach is challenging, as even then the optimization landscape is bimodal.Overall, the key difference in philosophy, setting and ultimately results, comes from our focus on optimization: our method is a study in the power of local (first and second-order) optimization methods.Question: Does the paper discuss the limitations of the work performed by the authors? Answer: [No] Justification: The paper focuses primarily on presenting a theoretical result. As such, it does not delve into empirical validations or practical implementations where limitations would typically be more relevant. The nature of the result is more abstract and formal, and thus, potential limitations about practical application, empirical robustness, or computational scalability are not directly addressed within the scope of this work. Guidelines:• The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper. • The authors are encouraged to create a separate "Limitations" section in their paper.• The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be. • The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with a few runs. In general, empirical results often depend on implicit assumptions, which should be articulated. • The authors should reflect on the factors that influence the performance of the approach.For example, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon. • The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size. • If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness. • While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome to describe the particular way they provide for reproducibility.In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results. 5. Open access to data and code Question: Does the paper provide open access to the data and code, with sufficient instructions to faithfully reproduce the main experimental results, as described in supplemental material? Answer: [Yes] Justification: We have provided the code we ran for our simulations in the supplementary material.Guidelines:• The answer NA means that paper does not include experiments requiring code. Justification: Since this is a theoretical result, we do not foresee any societal impact.Guidelines:• 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). 
preliminaries and notation: 
 In this section, we introduce essential definitions and notation for our work. We consider without loss of generality unit norm data points x 1 , . . . , x n ∈ R d , which we aim to project into k dimensions while preserving their norms with distortion at most ε = O( log n/k). Specifically, we seek matrices that satisfy the Johnson-Lindenstrauss guarantee:Definition 1 (Johnson-Lindenstrauss guarantee). The Johnson-Lindenstrauss guarantee (JL guarantee) states that for given dataset x 1 , . . . , x n ∈ R d and target dimension k, the distortion for all points does not exceed O( log n/k).To achieve this we define a linear mapping f (x) = Ax, where A ∈ R k×d . The JL Lemma guarantees that there exists a random linear mapping that achieves this projection with high probability:Lemma 1 (Distributional Johnson-Lindenstrauss Lemma). For ε, δ ∈ (0, 1) and k = O(log(1/δ)/ε 2 ), there exists a probability distribution D over linear functions f : R d → R k such that for every x ∈ R d :Pr f ∼D ∥f (x)∥ 2 2 ∈ (1 -ε)∥x∥ 2 2 , (1 + ε)∥x∥ 2 2 ≥ 1 -δ.There has been significant research aimed at improving the construction of these random mappings.In contrast to traditional algorithms, our approach proposes learning the linear mapping directly from the data, leveraging the inherent structure to surpass worst-case performance.Next, we give essential definitions for our optimization framework.Definition 2. A function f : R d → R is defined to be L-smooth, if for all x, y ∈ R d it satisfies:∥∇f (x) -∇f (y)∥ 2 ≤ L∥x -y∥ 2 .A function is called K-Hessian Lipschitz if for all x, y ∈ R d :∥∇ 2 f (x) -∇ 2 f (y)∥ 2 ≤ K∥x -y∥ 2 .Below, we give the definition for approximate stationarity.Definition 3. (Approximate stationarity). For a K-Hessian Lipschitz function f (•), we say that a point x * is a ρ-second-order stationary point (ρ-SOSP) if:∥∇f (x * )∥ 2 ≤ ρ and λ min (∇ 2 f (x * )) ≥ -Kρ.Notation. For vectors u, v we use ⟨u, v⟩ to denote their inner product and ∥u∥ 2 to denote the L 2 norm. For matrix M ∈ R k×d , we denote the element of the i th row and j th column by µ i,j and we use ∥M ∥ F to denote the Frobenius norm. For matrix M ∈ R k×d and σ 2 ∈ R + we use N (M , σ 2 ) to denote an k × d random Gaussian matrix where each element a i,j ∼ N (µ i,j , σ 2 ). We use ∇f and ∇ 2 f to denote the gradient and Hessian operators, respectively. 
the main results: 
 This section contains the full statement of our main theorems, and proof outlines. We organize the flow of this section according to our conceptual outline given in the introduction. In most cases, we defer the full proofs to the appendix.Our goal is to find a matrix that satisfies the Johnson-Lindenstrauss guarantee as given in Lemma 1.Consider the natural objective function of maximum distortion:h(A) = max x1,...,xn ∥Ax∥ 2 2 -1 . (1)Step 1: The first step tells us what will not work. In particular, direct optimization over the space of matrices cannot work. Our first result shows that minimizing this maximum distortion objective via a first or second-order method, is a doomed approach. In particular, we show that there exist instances which are bad local minima.Theorem 1. For all k > 1, there exists a family of matrices A k×k+1 which are strict local minima for the objective function of Equation 1 reachable from the origin. The achieved distortion is Ω(1) over a set of O(k 2 ) points, while there exist matrices yielding distortion O( log k/k) → 0.The proof of this is constructive. We construct a dataset and then show that a set of matrices reachable from the origin have large constant distortion, yet these points are locally unimprovable. The full proof can be found in Appendix A.1 □Step 2: The key idea towards our final result is to perform an optimization over an extended space: the space of parameters of random Gaussian solution samplers. We first define an optimization objective over this space, and then prove properties about the resulting landscape over the space of samplers.A Gaussian solution sampler is defined by its mean and variance. We only consider the case where all entries have the same variance. Thus, our new parameter space consists of pairs (M , σ 2 ), where M is a projection matrix, here interpreted as the mean of a Gaussian distribution, and σ 2 is the variance parameter. Given (M , σ 2 ), the solution sampler defined is simply: A ∼ N (M , σ 2 ).We note that our new parameter space has just one additional parameter than the ambient setting.Step 2A: We next must extend the maximum distortion objective above, to the space of random solution samplers we have defined. Consider the objective f * , defined as follows:EQUATIONThus, f * (M , σ 2 ) is the probability that a matrix sampled according to the corresponding Gaussian distribution will have maximum distortion at least ε. When we take ε = O( log n/k), our objective function f * gives us the probability that a Gaussian solution sampler fails to produce a matrix that meets the JL guarantee. Hence, a good sampler is one that makes f * small.We note that in the proof of the JL lemma in [Indyk and Motwani, 1998 ], the authors show that a matrix with Gaussian entries satisfies the JL guarantee with high probability. Thus, in particular, we know that taking M = 0, where 0 is a k × d zero matrix, and σ 2 = 1, gives a low objective value for f * .In the context of these definitions, therefore, our goal is to find a matrix M such that f * ( M , 0) has a low objective value. To do this, we now define a related objective value, which thanks to a regularization term, will allow the optimization algorithm to push us towards lower variance solutions. The technical challenge is then to show that we can control any deterioration of the JL guarantee of these lower variance solutions.We define our final objective function starting from f * defined above. First, we simplify the objective by applying a standard union bound and write a relaxed objective that sums for every point the probability that the point will have a norm outside the required bounds after the linear transformation.EQUATIONFor an appropriately chosen value of ε = O( log n/k), we have that no constraint is violated with probability greater than 1/(3n) and f (0, 1) < 1/3. The function f serves as an upper bound on the probability of generating a matrix that does not have the JL guarantee, effectively acting as a proxy for "bad" events. Next, we add an appropriate regularization term that penalizes high variance points. Our overall objective is thus:EQUATION)At our initialization point M = 0 and σ 2 = 1, the value of the regularization is 1/2, thus: g(0, 1) < 1/3 + 1/2 < 5/6. This is crucial, as it implies that following any decreasing path in g leads to points with a likelihood of a bad event being less than 1. Consequently, this convergence toward a solution sampler maintains a positive (and O(1)) probability of achieving a projection matrix that satisfies the JL guarantee. The next step provides our algorithm. After that, we characterize its fixed points.Step 2B: Algorithm 1 is a second-order descent algorithm consisting of two simple steps: At a given point x t = (M t , σ t ), if the gradient is sufficiently large, we take a gradient step. If the gradient is small, we consider the Hessian; if the smallest eigenvalue is sufficiently negative, we take a step in that direction of negative curvature; otherwise the algorithm terminates by reporting the mean parameter M t (see Lemma 4 for discussion on this final point). To prove correctness of the algorithm, we must show that any ρ-SOSP will have sufficiently small variance. The proof of correctness is given in Step 2C, and a bound on its running time in Step 3. We can call this algorithm recursively, to find the best distortion, using a simple routine given in Algorithm 2.We note that in principle, many first-order methods can also be used, for example, Perturbed Gradient Descent (PGD) which has been shown to converge to second-order stationary points fast [Jin et al., 2017] . We use a deterministic algorithm in our analysis to enable a straightforward derandomization of the JL lemma through the optimization of Equation 4.Algorithm 1 Hessian Descent.Require: ∇g, ∇ 2 g, ν = 1 L , h = 3 √ ρ K , L, K, ρ, M init = 0, σ 2 init = 1 1: t ← 0 2: x t ← (M init , σ 2 init ) 3: while true do 4: if ∥∇g(x t )∥ > ρ then 5: x t+1 ← x t -ν • ∇g(x t ) 6: else if ∥∇g(x t )∥ 2 ≤ ρ and λ min (∇ 2 g(x t )) < - √ Kρ then 7:u 1 ← the eigenvector corresponding to λ min (∇ 2 g(x t ))8:x t+1 ← x t + hu 1 9: else 10:return x t [0] = M t 11:end if 12:t ← t + 1 13: end whileStep 2C: Since we initialize at a good solution sampler and we use a descent algorithm on our objective, we know that we can never move to a bad sampler. But that is not enough for us. For we recall that our goal is not to find a good randomized algorithm, but rather to find a good (deterministic) JL matrix, via optimization. Thus we must show that we do not get trapped in any points that have non-zero variance.We accomplish this in several lemmas. First in Lemma 2 we show that stationary points must have zero variance. We then refine this in Lemma 3 and show that being in a ρ-second-order stationary point requires the variance to be very small. We need this in order to show we can escape from any point with sufficiently large variance. Finally, in Theorem 2 we show that our second-order algorithm Algorithm 2 An algorithm to find optimal distortion.Require: ∇g, ∇ 2 g, L, K, ρ, x initial , ε grid 1: min ε ← ∞ 2: min value ← ∞ 3: for each ε in ε grid do 4: M ε ← HESSIAN DESCENT(∇g, ∇ 2 g, ν, h, L, K, ρ) 5: value ← max distortion of M ε 6:if value < min value then 7:min value ← value 8:min ε ← ε 9:end if 10: end for 11: return min ε , min value will not get stuck at any point with large variance, and that once we are at a solution sampler with small enough variance, the mean parameter itself will enjoy (deterministically) the JL guarantee.In the following lemma, we show that points with non-zero variance cannot be local minima. Specifically, we show that for any given mean matrix and variance, there exists a nearby mean matrix with reduced variance that improves the objective value.Lemma 2. Let M ∈ R k×d , and σ > 0. Then for any γ ∈ [0, σ], there exists M ′ ∈ R k×d such that:• ∥M -M ′ ∥ F ≤ 2γ kd log 3 √ kd γ , • g(M ′ , σ 2 -γ 2 ) ≤ g(M , σ 2 ) -γ 2 /6.Proof Sketch: The proof of the lemma essentially is a small derandomization step, where we show that by taking a sufficiently small variance-reducing step, even if we deteriorate the JL guarantee (i.e., the function f increases), the decrease in the regularizer outweighs this increase, thereby decreasing the overall value of the objective, thus showing we could not have been at a local minimum.More specifically, the proof proceeds as follows. We begin with a Gaussian matrix A ∼ N (M , σ 2 ) and use the additivity property to partition it:A = A γ + A ′ , where A γ ∼ N (0, γ 2 ) and A ′ ∼ N (M , σ 2 -γ 2 ), representing small additive noise and the remainder of A, respectively.The core idea is to derandomize A γ to achieve a decreased objective value. We extend the definition of a bad event by constraining A γ to take values only within a specific range R. Within R, there must exist a realization of A γ , denoted as α γ , which at worst, only slightly increases the probability of failure due to the truncation of the Gaussian distribution tails. By choosing this specific value α γ , we effectively derandomize A γ . As we show in the appendix, this can result in an only slightly increased probability of a bad event.We then define M ′ = α γ + M and show that the regularization term ensures an overall decrease in the objective function. The full proof can be found in Appendix A.2. □In Lemma 2, we established that any point with non-zero variance cannot be a local minimum, as there always exists a nearby point with lower objective value. The next lemma addresses whether a descent direction can be found at each step. We prove that this is indeed the case, specifically demonstrating that the ρ-second-order stationary points of the objective function in Equation 4 correspond to points with approximately zero variance.Lemma 3. Consider x 1 , . . . , x n ∈ R d . Given target dimension k choose ε = O( log n/k).The ρ-second-order stationary points of the objective function in Eq. 4 impliesσ 2 < poly(n, k, d) • ρ O(1) .Proof Sketch: We establish the lemma by examining the behavior of the variance σ 2 at points approaching ρ-second-order stationarity under the objective function defined in Equation 4. While σ 2 is large we employ a series of incremental reductions, and we show we can continue in this manner until σ 2 is reduced at least to the claimed level.Using Taylor's theorem and the Lipschitzness of ∇ 2 g, we prove that at any point (M , σ 2 ), either the gradient will be large and thus progress will be made using first-order methods, or that the minimum eigenvalue of the Hessian will be negative and thus we can follow that direction to make progress. We then show that convergence to ρ-second-order stationary points gives us the desired result. Controlling the effect of the Lipschitz constant is a main challenge. The full proof can be found in Appendix A.3. □Step 2D: Since our Algorithm 1 finds an approximate ρ-SOSP, we need an additional result that gives us a stopping criterion once the variance is small enough, and simply use the mean with controlled deterioration of the JL guarantee. That is, instead of sampling from A ∼ N (M , σ 2 ), we can directly use M instead. This is why in line 10 of Algorithm 1, we simply return the parameter M t .Lemma 4. Given n unit vectors in R d and a target dimension k, choose ε = O( log n/k) such that distribution A ∼ N (M , σ 2 ) satisfies the JL guarantee with distortion ε with probability 1/6. Then using matrix M instead of sampling from A retains the JL guarantee with a threshold increased by at most poly(σ, 1/k).Proof Sketch: By assumption, A ∼ N (M , σ 2 ) satisfies 1/k∥Ax∥ 2 2 ∈ (1 -ε, 1 + ε) with probability at least 1/6. Next, we decompose A as A = M + Z with Z ∼ N (0, σ 2 ), and invoking the JL lemma, we chooseε 0 such that 1/k∥Zx∥ 2 2 ∈ [σ 2 (1 -ε 0 ), σ 2 (1 + ε 0 )] with probability at least 6/7. This choice ensures that there exists a region in the space such that the JL guarantee holds simultaneously for both A and Z. Our goal is to establish a bound on the distortion when using M instead of sampling from A. We show that ∥M x∥ 2 2 can be bounded with terms involving only ∥Ax∥ 2 2 and ∥Zx∥ 2 2 and derive an upper and lower bound on the distortion that must hold deterministically. The full proof can be found in Appendix A.4. □Putting the above results together, we have our first main result that shows the correctness of Algorithm 1: it will not get trapped at any point with a large variance, and once it does finally arrive at a point of small enough variance, the mean parameter satisfies the JL property. Theorem 2. Given n unit vectors in R d and a target dimension k, consider the corresponding function g of Equation 4for any ε ≥ C log n/k where C is a sufficiently large constant. Any ρ-SOSP, that is a pair of parameters (M , σ 2 ), of g reachable from the origin satisfies σ 2 < poly(n, k, d)• ρ O(1)and goes to 0 as ρ → 0. Moreover, when ρ < 1/poly(n, k, d), M satisfies the JL guarantee having distortion at most O(ε).Proof. We choose the parameter C such that g(0, 1) < 5/6. Using Lemmas 2 and 3 we find that a ρ-SOSP of g is a pair of parameters (M , σ 2 ), with σ 2 < poly(n, k, d) • ρ O(1) and g(M , σ 2 ) < g(0, 1). This implies that drawing sample from A ∼ N (M , σ 2 ) satisfies the JL guarantee with distortion at most ε with probability 1/6.Then, from Lemma 4, using the matrix M satisfies the JL guarantee with distortion O(ϵ). □Step 3: The above result proves correctness of Algorithm 1, but is qualitative: we have not proved how many steps are required. Below we give a quantitative result proving that we can minimize our objective function efficiently and learn a deterministic JL embedding in polynomial time, while incurring a minor increase in the distortion. Though as mentioned, we see the main contribution of our work as centering on the optimization formulation, we note that this theorem constitutes a novel approach to derandomizing the Gaussian JL transformation.Theorem 3. Given n unit vectors in R d and a target dimension k, consider any ε ≥ C log n/k where C is a sufficiently large constant. Then, running Algorithm 1 to deterministically optimize the objective function g of Equation 4for poly(n, k, d) steps returns a matrix M that satisfies the JL guarantee with distortion at most O(ε).Proof. The first part of this theorem relies on two auxiliary results: Lemma 5 (Sufficient Descent for Gradient Descent). If ∥∇g(x t )∥ 2 > ρ and L the smoothness of g, then for ν = 1 L and x t+1 = x t -ν • ∇g(x t ), we have g(x t+1 ) ≤ g(x t ) -νρ 2 2 . Lemma 6 (Sufficient Descent for Negative Curvature Descent). If ∥∇g(x t )∥ 2 ≤ ρ and λ min (∇ 2 g(x t )) < -√ Kρ, and K the Hessian Lipschitz parameter, then for h = 3 √ ρ K and x t+1 = x t + hu 1 , where u 1 corresponds to the eigenvector of the minimum eigenvalue, we haveg(x t+1 ) ≤ g(x t ) -3ρ 1.5 4 √ K .The proofs of both lemmas are in Appendices A.5, A.6, respectively. Given these two results, now standard optimization techniques show that Algorithm 1 finds a ρ-SOSP in O(1/ρ 1.5 ) steps.According to Theorem 2, choosing ρ < 1/poly(n, k, d), implies that a ρ-SOSP for g is a pair of parameters (M , σ 2 ) with M satisfying the JL guarantee with distortion at most O(ε). Therefore, running Algorithm 1 for poly(n, k, d) steps returns a matrix M which satisfies the JL guarantee with distortion at most O(ε). □ 
simulations: 
 We empirically validate our theoretical results, demonstrating that first and second-order information suffices to learn high-quality embeddings. Our goal is to identify an optimal matrix M that produces embeddings satisfying the JL guarantee (see Definition 1) while minimizing distortion. We show that our method generates embeddings with significantly lower distortion compared to those obtained using the Gaussian construction of the JL lemma.We generate a unit norm dataset with n = 100 data points in d = 500 dimensions, and target dimension k = 30 dimensions. We aim to minimize the expected maximum distortionE A∼N (M ,σ 2 ) [h(A)](see Equation 1). We employ the first-order optimization algorithm Adam [Kingma and Ba, 2014] over 5000 iterations and compare our method against the average and minimum distortions over 1000 trials using the baseline matrix Z ∼ N (0, 1) (left plot of Figure 1 ). To calculate the distortions of our method, we sample from the updated mean matrix and variance at each iteration.As predicted by our theoretical analysis, we demonstrate that while the Gaussian randomized construction achieves satisfactory distortion levels, our method converges to a high-quality deterministic solution that nearly eliminates distortion (both plots of Figure 1 ). At the conclusion of the procedure, we calculate the distortion using the resultant mean matrix M . Interestingly, the results demonstrate that ∥M x∥ 2 2 ≈ 1 (i.e. almost 0 distortion), compared to approximately 1 and 0.6 for the average and minimum distortions of the random construction, respectively. We plot the progression of the distortions and variance over the iterations in the plots of Figure 1 . This evidence highlights the effectiveness of our methodology in practice, showcasing the advantage of integrating data structure into dimensionality reduction for more accurate embeddings with provable guarantees.Further details can be found in Appendix B. Figure 1: Plot of the distortion obtained through optimization over 5000 iterations vs the average distortion using a random Gaussian matrix (left plot), and the progression of variance over the same number of iterations (right plot). To calculate the distortions' progression with our method, we sample from the updated mean matrix and variance at each iteration and compute the distortion. We remark that the distortion plotted is a proxy for our objective in Equation 4. We observe that our optimization-based approach converges to a deterministic solution sampler. By using the mean matrix M , we achieve nearly optimal distortion, where |M x| ≈ |x|. 
conclusion: 
 In this work, we have demonstrated that optimization used directly in the sampler space can find a deterministic JL-quality embedding. While our initial focus has been on the relatively simple random construction of the Johnson-Lindenstrauss Lemma, there exist more complex randomized constructions that may offer greater flexibility. We believe that the methods we have introduced could find applicability far outside the JL setting. 
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: Since this is a theoretical result, we do not think it has any risk of misuse 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: [NA] Justification: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. 
new assets: 
 Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?Answer: [NA] Justification: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: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. • According to the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector. 
institutional review board (irb) approvals or equivalent for research with human subjects: 
 Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?Answer: [NA] Justification: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: A.1 Proof of Theorem 1Theorem. For all k > 1, there exists a family of matrices A k×k+1 which are strict local minima for the objective function of Equation 1 reachable from the origin. The achieved distortion is Ω(1) over a set of O(k 2 ) points, while there exist matrices yielding distortion O( log k/k) → 0.Proof. We show that trying to minimize the maximum distortion of points directly can lead to getting stuck in bad local minima. Specifically, there exists a set of points such that the family of matrices F ⊆ R k×(k+1) , which consist of unitary matrices scaled by 2 in the first k × k block and a zero vector in the last column, is a local minimum, and any matrix A ′ with ∥A ′ -A∥ 2 F < 1/poly(k) for some matrix A ∈ F results in worse maximum distortion.Consider any matrix A ′ and its closest matrix A ∈ F. A ′ can be written as[A ′ [k×k] |v] where A ′ [k×k]represents the first k × k block of A ′ and v its last column. If we compute the polar decomposition of A ′[k×k] = U P , where U is a unitary matrix and P is a positive semidefinite symmetric matrix, the closest matrix A to A ′ can be written as [2U |0]. By rotating the space with U -1 , we assume, without loss of generality, that U = I and thus A ′[k×k] is a positive definite symmetric matrix. Dataset Construction: Define vectors e i ∈ R k for i = 1, . . . , k. The first k components of our vectors will be either e i or e i + e j with i ̸ = j. For each such vector we have 4 different versions for the last coordinate, denoted by last(x), is given by: last(x) ∈ {+ √ 15∥x∥ 2 , - √ 15∥x∥ 2 + √ 7/3∥x∥ 2 , - √ 7/3∥x∥ 2 }Thus, the data points are constructed as x = (x, last(x)). Essentially, last(x) is the completion of the x depending on whether we want tight large distortion or tight small distortion.We analyze ∥Ax∥ 2 2 -∥A ′ x∥ 2 2 . Using the Taylor expansion:∥Ax∥ 2 2 -∥A ′ x∥ 2 2 = 2(Ax) T (A -A ′ )x + O(∥A ′ -A∥ 2 F ) Since ∥A -A ′ ∥ 2 F = ∥P -2I∥ 2 F + ∥v∥ 2 2 , let ∥A -A ′ ∥ 2 F = γ 2 .We consider two cases based on the magnitude of ∥v∥ 2 2 . Case 1: ∥v∥ 2 2 ≥ γ 2 /2 In this case, there exists at least one i such that v 2 i ≥ γ 2 /(2k). We consider:2(Ax) T (A -A ′ )x = x T (P -2I)x + x T v last(x)Choosing x = e i , we encounter two sub-cases, depending on the sign of x T (P -2I)x. In both sub-cases, the term 2(Ax) T (A -A ′ )x dominates the higher-order term. This holds because:|x T (P -2I)x + x T vlast(x)| > |x T (P -2I)x + √ 7/3γ| > γ 2 = ∥A -A ′ ∥ 2 F .Then, we have:-If x T (P -2I)x < 0, select the negative value of the tight large distortion case for last(x) = -√ 7/3∥e i ∥ to ensure ∥Ax∥ 2 2 < ∥A ′ x∥ 2 2 , resulting in worse distortion. That is, we expected norm 4/3, but instead using matrix A, we got norm 2, and any A ′ nearby will only increase gap.-If x T (P -2I)x ≥ 0, select the positive value of the tight small distortion case for last(x) = √ 15∥e i ∥ 2 to ensure ∥Ax∥ 2 > ∥A ′ x∥ 2 , again resulting in worse distortion. resulting in worse distortion. In this case, we expected norm 4, but instead we got norm 2, and any A ′ nearby will only increase gap.Case 2: ∥v∥ 2 2 < γ 2 /2 This implies that ∥P -2I∥ 2 F ≥ γ 2 /2. Again, we encounter two sub-cases. Either an entry on the diagonal of P -2I has mass ≥ γ/2k 2 or an off-diagonal entry does. We shall prove that in each case, the first-order term dominates.-Diagonal Entry: If |(P -2I) i,i | ≥ γ/2k 2 , choose x = e i . Then:|x T (P -2I)x + x T vlast(x)| > |γ/2k 2 + x T vlast(x)| > γ 2 = ∥A -A ′ ∥ 2 F .We note here that we always match the sign of last(x) to that of x T (P -2I)x to increase the error. Now, we can follow the same logic as in Case 1 to show that there always exists a point with worse distortion.-Off-Diagonal Entry: If |(P -2I) i,j | ≥ γ/2k 2 for some i ̸ = j, choose x = e i + e j . Since P is symmetric, we are adding the contribution of two elements, and even if the diagonal elements have the opposite sign, by definition of this sub-case, their magnitude does not suffice to make the norm of the first-order term negligible:|x T (P -2I)x + x T vlast(x)| > |γ/2k 2 + x T vlast(x)| > γ 2 = ∥A -A ′ ∥ 2 F .We can follow the same logic as in Case 1 to show that there always exists a point with worse distortion.In summary, in all cases, the first-order term x T (P -2I)x + x T vlast(x) dominates over higher-order terms O(∥A -A ′ ∥ 2 F ), and thus the sign is always determined by the first term. We have shown that there always exists a point x causing increased distortion, thus confirming that any perturbation ∥A ′ -A∥ 2 F < 1/poly(k) results in worse maximum distortion. Thus, this family of matrices is a local minimum with constant distortion. □A.2 Proof of Lemma 2Lemma. Let M ∈ R k×d , and σ > 0 and consider a random matrix A ∼ N (M , σ 2 ). Then for any γ ∈ [0, σ], there exists M ′ ∈ R k×d such that:• ∥M -M ′ ∥ F ≤ 2γ kd log 3 √ kd γ • g(M ′ , σ 2 -γ 2 ) ≤ g(M , σ 2 ) -γ 2 /6Proof. Consider a random matrix A ∼ N (M , σ 2 ). Let γ ∈ [0, σ), then using the additivity property of Gaussian distributions, we decomposeA into A γ ∼ N (0, γ 2 ) and A ′ ∼ N (M , σ 2 -γ 2 ), such that A = A γ + A ′ .We note that the case for γ 2 = σ 2 follows from the same analysis but takingA ′ = M deterministically.We extend the definition of a "bad event", and consider a failure when ∥Ax∥ 2 falls outside the range . This extension results in an increased probability of a bad event by γ 2 /3. Denote by a i,j the j th element of the i th row of matrix A γ , then:R 1 := k(1 -ε), k(1 + ε) forPr   |a i,j | ≥ 2γ log √ kd γ/3 , ∀i = 1, . . . , k, ∀i = 1, . . . , d   = k i=1 d j=1   |a 1,i | ≥ 2γ log √ kd γ/3   ≤ k i=1 d j=1 2 exp   - 4γ 2 log √ kd γ/3 2γ 2   = k i=1 d j=1 2 exp log γ 2 /9 kd = k i=1 d j=1 2 γ 2 /9 kd = 2γ 2 /9 < γ 2 /3.We note that the inequality holds due to standard Gaussian properties.We prove that there exists a realization for A γ ∼ N (0, γ 2 ) within the range R 2 , which leads to an improvement in the objective function. To see this, we express Equation 3 as follows:EQUATIONLet δ denote the current value of f (M , σ 2 ) and consider the following expression of the extended definition of a bad event:EQUATIONWe lower bound Equation 6 by the minimum of the function inside the expectation to show that there must exist a realization of A γ that lies in the range R 2 , and increases the original objective function from δ to δ + γ 2 3 :min A γ ∈R2 n j=1 Pr A ′ [∥A γ x j + A ′ x j ∥ 2 ̸ ∈ R 1 | A γ ] < E A γ   n j=1 Pr A ′ [∥A γ x j + A ′ x j ∥ 2 ̸ ∈ R 1 | A γ ] + 1[A γ ̸ ∈ R 2 ]   < δ + γ 2 /3. Denote α γ := arg min A γ ∈R2 n j=1 Pr A ′ [∥A γ x j + A ′ x j ∥ 2 ̸ ∈ R 1 | A γ ], then A = α γ + A ′ ∼ N (M ′ , σ 2 -γ 2 ), where M ′ = α γ + M . Consequently, this means that:f (M ′ , σ 2 -γ 2 ) = n j=1 Pr A [∥Ax j ∥ 2 ̸ ∈ R 1 ] < δ + γ 2 /3.However, since we are working with the regularized objective in Equation 4, reducing σ 2 to σ 2 -γ 2 means that we have an overall decrease of the objective. Proof. Assume that we have reached a point v = (M , σ 2 ). Our analysis proceeds in two distinct cases based on the magnitude of σ 2 . If the variance exceeds a specified threshold, we employ a series of incremental reductions. This procedure continues until the variance reaches a sufficiently low threshold, facilitating a transition to a scenario wherein the variance can be reduced to zero in a single step.g(M ′ , σ 2 -γ 2 ) -g(M , σ 2 ) < (δ + γ 2 /3 -γ 2 /2) -δ < γ 2 /3 -γ 2 /2 < -γ 2 /6. □ A.3 Proof of Lemma 3 Lemma. Consider x 1 , . . . , x n ∈ R d . Given target dimension k choose ε = O( log n/k).Case 1: σ 2 > 2 -5 (Kd 3/2 log(kdK) 3/2 ) -2 . Choose γ 2 = 2 -5 (Kd 3/2 log(kdK) 3/2 ) -2 . According to Lemma 2, there exists a neighboring point M ′ and we define v ′ = (M ′ , σ 2 -γ 2 ) which reduces the objective function. Then, the distance between v ′ and v is given by:∥v ′ -v∥ 2 = 4γ 2 d log √ kd γ/3 + γ 4 .Denote K the Hessian Lipschitz constant, then, from Lemma 2 and Taylor's theorem on the Lipschitzness of ∇ 2 g [Nesterov, 2013] we get:g(v ′ ) -g(v) -∇g(v) • (v ′ -v) - 1 2 (v ′ -v) • ∇ 2 g(v) • (v ′ -v) ≤ K 6 ∥v ′ -v∥ 3 2 ∇g(v) • (v ′ -v) + 1 2 (v ′ -v) • ∇ 2 g(v) • (v ′ -v) - K 6 ∥v ′ -v∥ 3 2 ≤ g(v ′ ) -g(v) ∇g(v) • (v ′ -v) + 1 2 (v ′ -v) • ∇ 2 g(v) • (v ′ -v) - K 6 ∥v ′ -v∥ 3 2 ≤ - γ 2 6 .This implies that either∇g(v) • (v ′ -v) < - γ 2 12 =⇒ ∥∇g(v)∥ 2 > γ 12 4d log √ kd γ/3 + γ 2 , or (v ′ -v) • ∇ 2 g(v) • (v ′ -v) - K 6 ∥v ′ -v∥ 3 2 < - γ 2 12 .The last inequality implies thatEQUATION)The last inequality follows from the choice of γ 2 we made at the beginning.This guarantees that in each step either the gradient will be large and thus progress will be made using first-order methods or that the minimum eigenvalue of the Hessian will be negative and thus there exists a direction which we can follow by a second-order method.Case 2: σ 2 ≤ 2 -5 (Kd 3/2 log(kdK) 3/2 ) -2 .Assume that after several iterations, we have reached a point v = (M , σ 2 ). Given that the inequality in 7 holds for all γ 2 ≤ σ 2 , we select the maximum allowable reduction. This allows us to reduce the variance to zero in a single step. Based on Lemma 2, there exists a neighboring point v ′ = (M ′ , 0) which reduces the objective function and by following the steps as in the previous case, we get:∥∇g(v)∥ 2 > σ 12 4d log √ kd σ/3 + σ 2 . or λ min (∇ 2 g(v)) < - 1 24 4d log √ kd σ/3 + σ 2 .This means that, convergence to a ρ-second-order stationary point for g implies:σ 12 4d log √ kd σ/3 + σ 2 < ρ, 1 24 4d log √ kd σ/3 + σ 2 < Kρ.We will use the second inequality. First, we observe that:1 24 12 √ kd σ + σ 2 ≤ 1 24 4d log √ kd σ/3 + σ 2Therefore, we have:1 √ Kρ < 288 √ kd σ + 24σ 2Because of our analysis we have that at all times during the optimization σ 2 < 1 (refer to the regularized objective function in Equation 4. In addition to that, we assume that d ≥ 2, k ≥ 1, which are both natural assumptions. This in turn implies that:1 √ Kρ < 576 √ kd σ 1 576d √ Kρ < 3 √ kd σ 1 576d √ Kρkd < 1 σ σ 2 < 576 2 Kρkd 3 .Finally, substituting in for the Lipschitz constant K, we get:σ 2 < poly(n, k, d) • ρ O(1) .This implies that σ 2 can become arbitrarily small for an appropriate choice of ρ.A.4 Proof of Lemma 4Lemma. Given n unit vectors in R d and a target dimension k, choose ε = O( log n/k) such that distribution A ∼ N (M , σ 2 ) satisfies the JL guarantee with distortion ε with probability 1/6. Then using matrix M instead of sampling from A retains the JL guarantee with a threshold increased by at most poly(σ, 1/k).Proof. We start with the assumption that 1 k ∥Ax∥ 2 2 ∈ (1 -ε, 1 + ε) with probability at least 1 6 . Expressing A as A = M + Z where Z ∼ N (0, σ 2 ), and from the JL lemma, we can select ε 0 such that1 k ∥Zx∥ 2 2 ∈ [σ 2 (1 -ε 0 ), σ 2 (1 + ε 0 )],with probability at least 6 7 . This ensures there exists an overlap where both inequalities for A and Z hold simultaneously.Our goal is to determine how much worse the distortion becomes when using M instead of sampling from the distribution A.Using the triangle inequality we have:1 k ∥M x∥ 2 = 1 k ∥M x + Zx -Zx∥ 2 ≤ 1 k ∥M x + Zx∥ 2 + 1 k ∥Zx∥ 2 = 1 k ∥Ax∥ 2 + 1 k ∥Zx∥ 2 ,which by squaring both sides and using the JL guarantee for A and Z, we obtain:1 k ∥M x∥ 2 2 ≤ 1 k ∥Ax∥ 2 2 + 2 k 2 ∥Ax∥ 2 ∥Zx∥ 2 + 1 k ∥Zx∥ 2 2 ≤ 1 + ε + 2σ k √ 1 + ε √ 1 + ε 0 + σ 2 (1 + ε 0 ) ≤ 1 + ε + 2 √ 2σ k √ 1 + ε + 2σ 2 .For the lower bound, using the Cauchy-Schwarz inequality and the JL guarantee for A and Z, we have:1 k ∥M x∥ 2 2 ≥ 1 2k ∥M x + Zx∥ 2 2 - 1 k ∥Zx∥ 2 2 = 1 2k ∥Ax∥ 2 2 - 1 k ∥Zx∥ 2 2 ≥ 1/2(1 -ε) -σ 2 (1 + ε 0 ) ≥ 1/2(1 -ε) -σ 2 ,Combining these results, we observe that replacing A with M maintains the JL guarantee with an increased distortion threshold, bounded by at most poly(σ, 1/k).A.5 Proof of Lemma 5Proof. We use Taylor's theorem on the smoothness of g and get:Lemma. If ∥∇g(x t )∥ 2 > ρ, then for ν = 1 L and x t+1 = x t -ν • ∇g(x t ), we have g(x t+1 ) ≤ g(x t ) -νρ 2 2 . g(x t+1 ) ≤ g(x t ) + ⟨∇g(x t ), x t+1 -x t ⟩ + L 2 ∥x t+1 -x t ∥ 2 2 = g(x t ) -ν∥∇g(x t )∥ 2 2 + Lν 2 2 ∥∇g(x t )∥ 2 2 = g(x t ) -1 - 1 2 Lν νρ 2 .We can use ν = 1/L to get the following:g(x t+1 ) ≤ g(x t ) - 1 2 νρ 2 . □ A.6 Proof of Lemma 6 Lemma. If ∥∇g(x t )∥ 2 ≤ ρ and λ min (∇ 2 g(x t )) < - √ Kρ, then for h = 3 √ ρK and x t+1 = x t + hu 1 , where u 1 corresponds to the eigenvector of the minimum eigenvalue, we have g(x t+1 ) ≤ g(x t ) -3ρ 1.5 4 √ K .Proof. We use Taylor's theorem on the Lipschitzness of ∇ 2 g [Nesterov, 2013] and get:g(x t+1 ) ≤ g(x t ) + ⟨∇g(x t ), x t+1 -x t ⟩ + 1 2 ⟨∇ 2 g(x t )(x t+1 -x t ), (x t+1 -x t )⟩ + K 6 ∥x t+1 -x t ∥ 3 = g(x t ) + h⟨∇g(x t ), u 1 ⟩ + h 2 λ 1 2 + h 3 K 6 ≤ g(x t ) + hρ -h 2 Kρ + h 3 K 6 .We can use h = 3 √ ρ √ K to get the following:g(x t+1 ) ≤ g(x t ) - 3ρ √ ρ 2 √ K .We explore the behavior of a distortion optimization process using by minimizing the expected maximum distortion in a given dataset through a series of optimization descent steps. We show that while the Gaussian randomized construction can achieve good enough distortion, our method goes beyond that and achieves almost zero distortion, by taking into account the structure of the data.B.1 Proxy for the Objective Function.To do this, we use the maximum distortion (Equation 1) with a proxy of our objective function (Equation 4) and we aim to minimize:EQUATION)To do this we use the gradient of Equation 8 with respect to the parameters of interest θ = (M , σ 2 ):EQUATION)We can approximate the expectation by taking y 1 , . . . , y N samples drawn from N (M , σ 2 ) and using Monte Carlo sampling we get the approximate gradient:∇ (M ,σ 2 ) f (M , σ 2 ) ≈ 1 N N i=1 h(y i ).B.2 Methodology of the Simulation.We generate a unit norm synthetic dataset of n = 100 data points in d = 500 dimensions and our goal is to project these into k = 30 dimensions while minimizing f . We run our optimization for 5000 iterations, using the Adam optimizer [Kingma and Ba, 2014] with a learning rate of 0.01 and batch size N = 20. At every iteration, we calculate the maximum distortion and store it. We demonstrate that, through this procedure, the model consistently reduces the distortion and the variance, thus converging to a deterministic solution sampler.We compare our method with the Gaussian random construction, that is we draw Z ∼ N (0, 1) and calculate ∥Zx∥ 2 2 . To have a fair evaluation we draw 1000 such matrices and calculate the mean and the minimum distortion obtained.Our results show that we consistently learn a matrix with close to optimal distortion, that is ∥M x∥ 2 2 ≈ 1, while the randomized construction achieves an average value of ∥Z avg x∥ 2 2 ≈ 2 and minimum value ∥Z min x∥ 2 ≈ 1.6.Denote by µ 1 , . . . , µ k , the rows of matrix M and A 1 , . . . , A k the rows of the Gaussian random matrix A. Then we can write the objective function:EQUATION)where f k,δ J , F k,δj are the pdf and cdf of χ 2 k , the non-central chi-squared distribution with k degrees of freedom andδ j = ⟨µ 1 , x j ⟩ 2 + • • • + ⟨µ k , x j ⟩ 2 σ 2as the non centrality parameter, respectively.Note that, instead of considering the k × d mean variables directly, it is simpler to reframe the problem. Specifically, we can view the function g in terms of the inner product variables v 1 , . . . , v k and σ 2 , where v i = ⟨µ i , x⟩ and µ i is the i-th row of the matrix M .Therefore, our problem reduces to proving the Lipschitz continuity of the Gradient and Hessian continuity with respect to these new variables. This approach simplifies the calculations significantly.To establish that the Gradient and Hessian is Lipschitz continuous, we examine the case for j = 1.The extension to j = n can be handled through summation. Let τ = σ 2 . Then, we have δ= ∥v∥ 2 τ = v 2 1 + • • • + v 2 k τ. Our function from Equation 10becomesg(v 1 , . . . , v k , τ ) = 1 + F k,δ k(1 -ε) τ -F k,δ k(1 + ε) τ + τ /2, where F k,δ (x) = e -δ/2 ∞ j=0 (δ/2) j j! Q(x; k + 2j), with Q(x; k) = γ(k/2,x/2) Γ(k/2)and γ(y, t) is the lower incomplete gamma function.Notice that taking the derivative gives an extra part with increased degrees of freedom∂F δ,k (x) ∂δ = -1/2F δ,k (x) + 1/2F δ,k+2 (x).C.1 Gradient Lipschitzness.Denote by D the elements of ∇g, and note that in this section to simplify notation we use ∥ • ∥ to represent ∥ • ∥ 2 . Then the derivatives are:D vi (v 1 , . . . , v k , τ ) = - v i τ F k,δ k(1 -ε) τ + v i τ F k+2,δ k(1 -ε) τ + v i τ F k,δ k(1 + ε) τ - v i τ F k+2,δ k(1 + ε) τ . D τ (v 1 , . . . , v k , τ ) = ∥v∥ 2 2τ 2 F k,δ k(1 -ε) τ - ∥v∥ 2 2τ 2 F k+2,δ k(1 -ε) τ - k(1 -ε) τ 2 f k,δ k(1 -ε) τ - ∥v∥ 2 2τ 2 F k,δ k(1 + ε) τ + ∥v∥ 2 2τ 2 F k+2,δ k(1 + ε) τ + k(1 + ε) τ 2 f k,δ k(1 + ε) τ .To prove Gradient Lipschitzness of g we can bound the Frobenius norm of the Hessian ∇ 2 g. Denote by D 2 the elements of ∇ 2 g. Below we calculate and bound all the derivatives:D 2 vi,vi g(v 1 , . . . , v k , τ ) = - 1 τ + v 2 i τ 2 F k,δ k(1 -ε) τ + - 2v 2 i τ 2 + 1 τ F k+2,δ k(1 -ε) τ + v 2 i τ 2 F k+4,δ k(1 -ε) τ + 1 τ - v 2 i τ 2 F k,δ k(1 + ε) τ + 2v 2 i τ 2 - 1 τ F k+2,δ k(1 + ε) τ - v 2 i τ 2 F k+4,δ k(1 + ε) τ .Here we used the triangle inequality and the fact that the cdf is bounded by 1. Notice how there is no dependency on ∥v∥ since the range we are integrating over with the cumulative distribution functions is independent of ∥v∥. Consequently, if ∥v∥ is large, the probability becomes exponentially small. Thus, we can bound this by:|D 2 vi,vi g(v 1 , . . . , v k , τ )| ≤ 4 τ + 8 τ 2 ≤ O 1 τ 2 .Next, we have:D 2 vi,vj g(v 1 , . . . , v k , τ ) = v i v j τ 2 F k,δ k(1 -ε) τ - 2v i v j τ 2 F k+2,δ k(1 -ε) τ + v i v j τ 2 F k+4,δ k(1 -ε) τ - v i v j τ 2 F k,δ k(1 + ε) τ + 2v i v j τ 2 F k+2,δ k(1 + ε) τ - v i v j τ 2 F k+4,δ k(1 + ε) τ .Thus, we can bound this by:|D 2 vi,vj g(v 1 , . . . , v k , τ )| ≤ 8v i v j τ 2 ≤ O 1 τ 2 .Next, we have:D 2 vi,τ g(v 1 , . . . , v k , τ ) = v i τ 2 - v i ∥v∥ 2 τ 3 F k,δ k(1 -ε) τ + 2v i ∥v∥ 2 τ 3 - v i τ 2 F k+2,δ k(1 -ε) τ - v i ∥v∥ 2 τ 3 F k+4,δ k(1 -ε) τ + v i k(1 -ε) τ 3 f k,δ k(1 -ε) τ - v i k(1 -ε) τ 3 f k+2,δ k(1 -ε) τ - v i τ 2 + v i ∥v∥ 2 τ 3 F k,δ k(1 + ε) τ + - 2v i ∥v∥ 2 τ 3 + 2v i τ 2 F k+2,δ k(1 + ε) τ + v i ∥v∥ 2 τ 3 F k+4,δ k(1 + ε) τ - v i k(1 + ε) τ 3 f k,δ k(1 + ε) τ + v i k(1 + ε) τ 2 f k+2,δ k(1 + ε) τ .Thus, we can bound this by:|D 2 vi,τ g(v 1 , . . . , v k , τ )| ≤ 4v i τ 2 + 8v i ∥v∥ 2 τ 3 + 2v i k(1 -ε) τ 3 + 2v i k(1 + ε) τ 3 ≤ O k τ 3 .Next, we have:D 2 τ,τ g(v 1 , . . . , v k , τ ) = - ∥v∥ 2 τ 3 + ∥v∥ 4 4τ 4 F k,δ k(1 -ε) τ + - ∥v∥ 4 2τ 4 + ∥v∥ 2 τ 3 F k+2,δ k(1 -ε) τ + ∥v∥ 4 4τ 4 F k+4,δ k(1 -ε) τ + - ∥v∥ 2 k(1 -ε) 2τ 4 + 2k(1 -ε) τ 3 - (k(1 -ε)) 2 2τ 4 f k,δ k(1 -ε) τ + ∥v∥ 2 k(1 -ε) 2τ 4 f k+2,δ k(1 -ε) τ + k(1 -ε) τ 3 e -δ/2 ∞ j=0 (δ/2) j j! ((l + 2j)/2 -1)f k+2j k(1 -ε) τ + ∥v∥ 2 τ 3 - ∥v∥ 4 4τ 4 F k,δ k(1 + ε) τ + ∥v∥ 4 2τ 4 - ∥v∥ 2 τ 3 F k+2,δ k(1 + ε) τ - ∥v∥ 4 4τ 4 F k+4,δ k(1 + ε) τ + ∥v∥ 2 k(1 + ε) 2τ 4 - 2k(1 + ε) τ 3 + (k(1 + ε)) 2 2τ 4 f k,δ k(1 + ε) τ - ∥v∥ 2 k(1 + ε) 2τ 4 f k+2,δ k(1 + ε) τ - k(1 + ε) τ 3 e -δ/2 ∞ j=0 (δ/2) j j! ((l + 2j)/2 -1)f k+2j k(1 + ε) τ + 1 2 .Thus, we can bound this by:|D 2 τ,τ g(v 1 , . . . , v k , τ )| ≤ 4∥v∥ 2 τ 3 + 2∥v∥ 4 τ 4 + ∥v∥ 2 k(1 -ε) τ 4 + 3k(1 -ε) τ 3 + (k(1 -ε)) 2 2τ 4 + ∥v∥ 2 k(1 + ε) τ 4 + 3k(1 + ε) τ 3 + (k(1 + ε)) 2 2τ 4 ≤ O kε 2 τ 4 .Overall, for D 2 vi,vi g, we have:ξ vi,vi = max v1,...,v k ,τ ∥D 2 vi,vi g(v 1 , . . . , v k , τ )∥ ≤ min τ (1/τ 2 ).Similarly, for D 2 vi,vj g, we obtain ξ vi,vj ≤ min τ (1/τ 2 ) for D 2 vi,τ g, we obtain ξ vi,τ ≤ min τ (k/τ 3 ) and for D 2 τ,τ , we obtain ξ τ,τ ≤ min τ (k 2 ε 2 /τ 4 ). Substituting M and σ 2 back, and assuming that the minimum value for σ 2 we allow is σ 2 0 we have that the Frobenius norm of the Hessian is bounded from:∥∇ 2 g∥ 2 F ≤ k i=1 ξ 2 vi,vi + ξ 2 σ 2 ,σ 2 + 2 i̸ =j ξ 2 vi,vj + 2 k i=1 ξ 2 vi,σ 2 ≤ kξ 2 vi,vi + ξ 2 σ 2 ,σ 2 + 2(k × d -2k -1)ξ 2 vi,vj + 2kξ 2 vi,σ 2 .Therefore we get:∥∇ 2 g∥ F ≤ O k σ 4 0 + O k 2 ε 2 σ 8 0 + O dk σ 4 0 + O k σ 6 0 ≤ poly (k, ε, d, 1/σ 0 ) .Finally, since we are summing over n data points the smoothness Lipschitz constant is poly(n, k, ε, d, 1/σ 0 )C.2 Hessian Lipschitzness.To prove Hessian Lipschitzness of g we will use the definition and the third-order derivatives. Denote by D 3 the elements of ∇ 3 g. Below we calculate and bound all the third-order derivatives:D 3 vi,vi,vi g(v 1 , . . . , v k , τ ) = 3v i τ 2 - v 3 i τ 3 F k,δ k(1 -ε) τ + - 6v i τ 2 + 3v 3 i τ 3 F k+2,δ k(1 -ε) τ + 3v i τ 2 - 3v 3 i τ 3 F k+4,δ k(1 -ε) τ + v 3 i τ 3 F k+6,δ k(1 -ε) τ + - 3v i τ 2 + v 3 i τ 3 F k,δ k(1 + ε) τ + 6v i τ 2 - 3v 3 i τ 3 F k+2,δ k(1 + ε) τ + - 3v i τ 2 + 3v 3 i τ 3 F k+4,δ k(1 + ε) τ - v 3 i τ 3 F k+6,δ k(1 + ε) τ .Thus, we can bound this:|D 3 vi,vi,vi g(v 1 , . . . , v k , τ )| ≤ 24v i τ 2 + 16v 3 i τ 3 ≤ 1 τ 3 .Here we used the triangle inequality and the fact that the cdf is bounded by 1. Notice how there is no dependency on ∥v∥ since the range we are integrating over with the cumulative distribution functions is independent of ∥v∥. Consequently, if ∥v∥ is large, the probability becomes exponentially small.Next, we have:D 3 vi,vj ,vj g(v 1 , . . . , v k , τ ) = v i τ 2 - v i v 2 j τ 3 F k,δ k(1 -ε) τ + - 2v i τ 2 + 3v i v 2 j τ 3 F k+2,δ k(1 -ε) τ + v i τ 2 - 3v i v 2 j τ 3 F k+4,δ k(1 -ε) τ + v i v 2 j τ 3 F k+6,δ k(1 -ε) τ + - v i τ 2 + v i v 2 j τ 3 F k,δ k(1 + ε) τ + 2v i τ 2 - 3v i v 2 j τ 3 F k+2,δ k(1 + ε) τ + - v i τ 2 + 3v i v 2 j τ 3 F k+4,δ k(1 + ε) τ - v i v 2 j τ 3 F k+6,δ k(1 + ε) τ .Thus, we can bound this by:|D 3 vi,vj ,vj g(v 1 , . . . , v k , σ)| ≤ 8v i τ 2 + 16v i v 2 j τ 3 ≤ 1 τ 3 .Next, we have:D 3 vi,vj ,v k g(v 1 , . . . , v k , τ ) = - v i v j v k τ 3 F k,δ k(1 -ε) τ + 3v i v j v k τ 3 F k+2,δ k(1 -ε) τ - 3v i v j v k τ 3 F k+4,δ k(1 -ε) τ + v i v j v k τ 3 F k+6,δ k(1 -ε) τ + v i v j v k τ 3 F k,δ k(1 + ε) τ - 3v i v j v k τ 3 F k+2,δ k(1 + ε) τ + 3v i v j v k τ 3 F k+4,δ k(1 + ε) τ - v i v j v k τ 3 F k+6,δ k(1 + ε) τ .Thus, we can bound this by:|D 3 vi,vj ,v k g(v 1 , . . . , v k , τ )| ≤ 16v i v j v k τ 3 ≤ 1 τ 3 .Next, we have:D 3 vi,τ,vi g(v 1 , . . . , v k , τ ) = 1 τ 2 - ∥v∥ 2 τ 3 - 3v 2 i τ 3 + v 2 i ∥v∥ 2 τ 4 F k,δ k(1 -ε) τ + 6v 2 i τ 3 - 3v 2 i ∥v∥ 2 τ 4 + 2∥v∥ 2 τ 3 - 1 τ 2 F k+2,δ k(1 -ε) τ + 3v 2 i ∥v∥ 2 τ 4 - 3v 2 i τ 3 - ∥v∥ 2 τ 3 F k+4,δ k(1 -ε) τ - v 2 i ||v∥ 2 τ 4 F k+6,δ k(1 -ε) τ + k(1 -ε) τ 3 - v 2 i k(1 -ε) τ 4 f k,δ k(1 -ε) τ + - k(1 -ε) τ 3 + 2v 2 i k(1 -ε) τ 4 f k+2,δ k(1 -ε) τ - v 2 i k(1 -ε) τ 4 f k+4,δ k(1 -ε) τ + - 1 τ 2 + ∥v∥ 2 τ 3 + 3v 2 i τ 3 - v 2 i ∥v∥ 2 τ 4 F k,δ k(1 + ε) τ + - 6v 2 i τ 3 + 3v 2 i ∥v∥ 2 τ 4 - 2∥v∥ 2 τ 3 + 1 τ 2 F k+2,δ k(1 + ε) τ + - 3v 2 i ∥v∥ 2 τ 4 + 3v 2 i τ 3 + ∥v∥ 2 τ 3 F k+4,δ k(1 + ε) τ + v 2 i ||v∥ 2 τ 4 F k+6,δ k(1 + ε) τ + - k(1 + ε) τ 3 + v 2 i k(1 + ε) τ 4 f k,δ k(1 + ε) τ + k(1 + ε) τ 3 - 2v 2 i k(1 + ε) τ 4 f k+2,δ k(1 + ε) τ + v 2 i k(1 + ε) τ 4 f k+4,δ k(1 + ε) τ .Thus, we can bound this by:|D 3 vi,τ,vi g(v 1 , . . . , v k , τ )| ≤ 4 τ 2 + 8∥v∥ 2 τ 3 + 24v 2 i τ 3 + 16v 2 i ∥v∥ 2 τ 4 + + 2k(1 -ε) τ 3 + 4v 2 i k(1 -ε) τ 4 + 2k(1 + ε) τ 3 + 4v 2 i k(1 + ε) τ 4 ≤ k τ 4 .Next, we have:D 3 vi,τ,vj g(v 1 , . . . , v k , τ ) = - 3v i v j τ 3 + v i v j ∥v∥ 2 τ 4 F k,δ k(1 -ε) τ + 6v i v j τ 3 - 3v i v j ∥v∥ 2 τ 4 F k+2,δ k(1 -ε) τ + 3v i v j ∥v∥ 2 τ 4 - 3v i v j τ 3 F k+4,δ k(1 -ε) τ - v i v j ∥v∥ 2 τ 4 F k+6,δ k(1 -ε) τ + 2v i v j k(1 -ε) τ 4 f k,δ k(1 -ε) τ - 4v i v j k(1 -ε) τ 4 f k+2,δ k(1 -ε) τ + 2v i v j k(1 -ε) τ 4 f k+4,δ k(1 -ε) τ + 3v i v j τ 3 - v i v j ∥v∥ 2 τ 4 F k,δ k(1 + ε) τ + - 6v i v j τ 3 + 3v i v j ∥v∥ 2 τ 4 F k+2,δ k(1 + ε) τ + - 3v i v j ∥v∥ 2 τ 4 + 3v i v j τ 3 F k+4,δ k(1 + ε) τ + v i v j ∥v∥ 2 τ 4 F k+6,δ k(1 + ε) τ - 2v i v j k(1 + ε) τ 4 f k,δ k(1 + ε) τ + 4v i v j k(1 + ε) τ 4 f k+2,δ k(1 + ε) τ - 2v i v j k(1 + ε) τ 4 f k+4,δ k(1 + ε) τ .Thus, we can bound this by:|D 3 vi,τ,vj g(v 1 , . . . , v k , τ )| ≤ 24v i v j τ 3 + 16v i v j ∥v∥ 2 τ 4 + 16v i v j k(1 -ε) τ 4 + 16v i v j k(1 + ε) τ 4 ≤ k τ 4 .Next, we have:D 3 vi,τ,τ g(v 1 , . . . , v k , τ ) = - 2v i τ 3 + 2||v|| 2 v i τ 4 - v i ||v|| 4τ 5 F k,δ k(1 -ε) τ + - 6||v|| 2 v i τ 4 + 3v i ∥v∥ 4 4τ 5 + 2v i τ 3 F k+2,δ k(1 -ε) τ + - v i ||v∥ 4 2τ 5 + 2v i ∥v∥ 2 τ 4 - v i ∥v∥ 4 4τ 4 F k+4,δ k(1 -ε) τ + v i ∥v∥ 4 4τ 4 F k+6,δ k(1 -ε) τ + - 3v i k(1 -ε) τ 4 + v i ∥v∥ 2 k(1 -ε) 2τ 5 + v i (k(1 -ε)) 2 τ 5 f k,δ k(1 -ε) τ + - v i ∥v∥ 2 k(1 -ε) τ 5 + 3v i k(1 -ε) τ 4 - v i (k(1 -ε)) 2 2τ 5 f k+2,δ k(1 -ε) σ 2 + v i ∥v∥ 2 k(1 -ε) 2τ 5 f k+4,δ k(1 -ε) τ - v i k(1 -ε) τ 4 e -δ/2 ∞ j=0 (δ/2) j j! ((l + 2j)/2 -1)f k+2j k(1 -ε) τ + v i k(1 -ε) τ 4 e -δ/2 ∞ j=0 (δ/2) j j! ((l + 2j + 2)/2 -1)f k+2j+2 k(1 -ε) τ + 2v i τ 3 - 2||v|| 2 v i τ 4 + v i ||v|| 4τ 5 F k,δ k(1 + ε) τ + 6||v|| 2 v i τ 4 - 3v i ∥v∥ 4 4τ 5 - 2v i τ 3 F k+2,δ k(1 + ε) τ + v i ||v∥ 4 2τ 5 - 2v i ∥v∥ 2 τ 4 + v i ∥v∥ 4 4τ 4 F k+4,δ k(1 + ε) τ - v i ∥v∥ 4 4τ 4 F k+6,δ k(1 + ε) τ + 3v i k(1 + ε) τ 4 - v i ∥v∥ 2 k(1 + ε) 2τ 5 - v i (k(1 + ε)) 2 τ 5 f k,δ k(1 + ε) τ + v i ∥v∥ 2 k(1 + ε) τ 5 - 3v i k(1 + ε) τ 4 + v i (k(1 + ε)) 2 2τ 5 f k+2,δ k(1 + ε) σ 2 - v i ∥v∥ 2 k(1 + ε) 2τ 5 f k+4,δ k(1 + ε) τ + v i k(1 + ε) τ 4 e -δ/2 ∞ j=0 (δ/2) j j! ((l + 2j)/2 -1)f k+2j k(1 + ε) τ - v i k(1 + ε) τ 4 e -δ/2 ∞ j=0 (δ/2) j j! ((l + 2j + 2)/2 -1)f k+2j+2 k(1 + ε) τ .Thus, we can bound this by:|D 3 vi,τ,τ g(v 1 , . . . , v k , τ )| ≤ 8v i τ 3 + 20||v|| 2 v i τ 4 + v i ||v|| 2τ 5 + 14v i ||v|| 4 4τ 4 + 8v i k(1 -ε) τ 4 + 2v i ∥v∥ 2 k(1 -ε) τ 5 + 3v i (k(1 -ε)) 2 2τ 5 + 8v i k(1 + ε) τ 4 + 2v i ∥v∥ 2 k(1 + ε) τ 5 + 3v i (k(1 + ε)) 2 2τ 5 ≤ k 2 ε 2 τ 5 .Next, we have:D 3 τ,τ,τ g(v 1 , . . . , v k , τ ) = 3∥v∥ 2 τ 4 - 3∥v∥ 4 2τ 5 + ∥v∥ 6 8τ 6 F k,δ k(1 -ε) τ + 4∥v∥ 4 τ 5 - 3∥v∥ 6 8τ 6 - 3∥v∥ 2 τ 4 F k+2,δ k(1 -ε) τ + 3∥v∥ 6 8τ 6 - 3∥v∥ 4 2τ 5 F k+4,δ k(1 -ε) τ - ∥v∥ 6 8τ 6 F k+6,δ k(1 -ε) τ + 4∥v∥ 2 k(1 -ε) τ 5 - ∥v∥ 4 k(1 -ε) 2τ 6 - 6k(1 -ε) τ 4 + 2(k(1 -ε)) 2 τ 5 - ∥v∥ 2 (k(1 -ε)) 2 4τ 6 f k,δ k(1 -ε) τ + ∥v∥ 4 k(1 -ε) τ 6 - 4∥v∥ 2 k(1 -ε) τ 5 + ∥v∥ 2 (k(1 -ε)) 2 4τ 6 f k+2,δ k(1 -ε) τ - ∥v∥ 4 k(1 -ε) 4τ 6 f k+4,δ k(1 -ε) τ + ∥v∥ 2 k(1 -ε) 2τ 5 - 3k(1 -ε) τ 4 + (k(1 -ε)) 2 2τ 5 e -δ/2 ∞ j=0 (δ/2) j j! (k + 2j -2)f k+2j k(1 -ε) τ - ∥v∥ 2 k(1 -ε) 2τ 5 e -δ/2 ∞ j=0 (δ/2) j j! (l/2 + j)f k+2+2j k(1 -ε) τ - k(1 -ε) τ 4 e -δ/2 ∞ j=0 (δ/2) j j! ((l + 2j)/2 -1) 2 f k+2j k(1 -ε) τ + - 3∥v∥ 2 τ 4 + 3∥v∥ 4 2τ 5 - ∥v∥ 6 8τ 6 F k,δ k(1 + ε) τ + - 4∥v∥ 4 τ 5 + 3∥v∥ 6 8τ 6 + 3∥v∥ 2 τ 4 F k+2,δ k(1 + ε) τ + - 3∥v∥ 6 8τ 5 + 3∥v∥ 4 2τ 5 F k+4,δ k(1 + ε) τ + ∥v∥ 6 8τ 6 F k+6,δ k(1 + ε) τ + - 4∥v∥ 2 k(1 + ε) τ 5 + ∥v∥ 4 k(1 + ε) 2τ 6 + 6k(1 + ε) τ 4 - 2(k(1 + ε)) 2 τ 5 + ∥v∥ 2 (k(1 + ε)) 2 4τ 6 f k,δ k(1 + ε) τ + - ∥v∥ 4 k(1 + ε) τ 6 + 4∥v∥ 2 k(1 + ε) τ 5 - ∥v∥ 2 (k(1 + ε)) 2 4τ 6 f k+2,δ k(1 + ε) τ + ∥v∥ 4 k(1 + ε) 4τ 6 f k+4,δ k(1 + ε) τ + - ∥v∥ 2 k(1 -ε) 2τ 5 + 3k(1 + ε) τ 4 - (k(1 + ε)) 2 2τ 5 e -δ/2 ∞ j=0 (δ/2) j j! (k + 2j -2)f k+2j k(1 + ε) τ + ∥v∥ 2 k(1 + ε) 2τ 5 e -δ/2 ∞ j=0 (δ/2) j j! (l/2 + j)f k+2+2j k(1 + ε) τ + k(1 + ε) τ 4 e -δ/2 ∞ j=0 (δ/2) j j! ((l + 2j)/2 -1) 2 f k+2j k(1 + ε) τ .Thus, we can bound this by:|D 3 τ,τ,τ g(v 1 , . . . , v k , τ )| ≤ 12∥v∥ 2 τ 4 + 7∥v∥ 4 τ 5 + ∥v∥ 6 τ 5 + 8∥v∥ 2 k(1 -ε) τ 5 + 7∥v∥ 4 k(1 -ε) 4τ 6 + 10k(1 -ε) τ 4 + 5(k(1 -ε)) 2 2τ 5 + ∥v∥ 2 (k(1 -ε)) 2 2τ 6 + 7∥v∥ 4 k(1 -ε) 4τ 6 + 8∥v∥ 2 k(1 + ε) τ 5 + 7∥v∥ 4 k(1 + ε) 4τ 6 + 10k(1 + ε) τ 4 + 5(k(1 + ε)) 2 2τ 5 + ∥v∥ 2 (k(1 + ε)) 2 2τ 6 + 7∥v∥ 4 k(1 + ε) 4τ 6 ≤ k 2 ε 2 τ 6 .Overall, for D 2 vi,vi g, we have:ρ vi,vi = max v1,...,v k ,τ ∥∇D 2 vi,vi g(v 1 , . . . , v k , τ )∥ ≤ max v1,...,v k ,τ {|D 3 vi,vi,vi g(v 1 , . . . , v k , τ )| + |D 3 vi,vi,vj g(v 1 , . . . , v k , τ )| + |D 3 vi,vi,τ g(v 1 , . . . , v k , τ )|} ≤ min τ (k/τ 4 ).Similarly, for D 2 vi,vj g, we obtain ρ vi,vj ≤ min τ (kε/τ 4 ) for D 2 vi,τ g, we obtain ρ vi,τ ≤ min τ (k 2 ε 2 /τ 5 ) and for D 2 τ,τ , we obtain ρ τ,τ ≤ min τ (k 2 ε 2 /τ 6 ). Substituting M and σ 2 back, and assuming that the minimum value for σ 2 we allow is σ 2 0 we have:∥H(M 1 , σ 2 1 ) -H(M 2 , σ 2 2 )∥ 2 = k i=1 g vi,vi (M 1 , σ 2 1 ) -g vi,vi (M 2 , σ 2 2 ) 2 + g σ 2 ,σ 2 (M 1 , σ 2 1 ) -g σ 2 ,σ 2 (M 2 , σ 2 2 ) 2 + i̸ =j 2 g vi,vj (M 1 , σ 2 1 ) -g vi,vj (M 2 , σ 2 2 ) 2 + 2 k i=1 g vi,σ 2 (M 1 , σ 2 1 ) -g vi,σ 2 (M 2 , σ 2 2 ) | 2 ≤ kρ 2 v1,v1 + ρ 2 σ 2 ,σ 2 + 2(k × d -2k -1)ρ 2 v1,v2 + 2kρ 2 v1,σ 2 (M 1 -M 2 ) 2 + (σ 2 1 -σ 2 2 ) 2 .Therefore, we get:∥H(M 1 , σ 2 1 ) -H(M 2 , σ 2 2 )∥ ≤ k 2 σ 8 0 + k 2 ε 2 σ 12 0 + dk 2 ε σ 8 0 + k 3 ε 2 σ 10 0 ∥(M 1 , σ 2 1 ) -(M 2 , σ 2 2 )∥ ≡ poly k, ε, d, 1 σ 0 ∥(M 1 , σ 2 1 ) -(M 2 , σ 2 2 )∥.Finally, since we are summing over n data points, the Hessian Lipschitz constant, is poly(n, k, ε, d, 1/σ 0 ).1. Claims Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope? Answer: [Yes] Justification: The abstract summarizes the findings presented in the paper, while the introduction summarizes the motivations behind our work and reviews prior research in the field.Guidelines:• The answer NA means that the abstract and introduction do not include the claims made in the paper. • The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers. • The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings. • It is fine to include aspirational goals as motivation as long as it is clear that these goals are not attained by the paper.38th Conference on Neural Information Processing Systems (NeurIPS 2024). 


Decision: This paper proposes a data-driven optimization approach for embedding that is in stark contrast to randomized methods such as the Johnson-Lindenstraus (JL). While solving the distance-preserving objective is intractable over the projection matrices, this paper considers an optimization program that searches a solution on a larger space and shows that the algorithm converges in polynomial steps.

The idea that optimizing over an expanded and carefully constructed space is novel. The results are sound, in the sense that the solution returned by the optimization program satisfies the JL guarantee; this is the first result of its kind which does not rely on random projection matrix. Empirically, it appears that the optimization-based approach enjoys lower distortion. That said, there are also weaknesses regarding the computational cost and empirical viability of the method as it relies on second order optimization.

Overall, the techniques presented in this paper may inspire a new line of research and represent a nice addition to NeurIPS program.

Reviews: review: paper_summary: paper_summary: This paper investigates the problem of using optimization-based approaches to learn Johnson Lindenstrauss(JL) embedding. The authors proposed a new framework to achieve the JL guarantee via optimization, instead of the traditional randomized methods. Similar with diffusion models, the authors proposed a novel method that uses optimization in an extended space of random Gaussian solution samplers, which circumvents direct optimization in non-convex landscape. Their approach uses second-order descent, gradually reduces the variance without increasing the expected distortion of the sampler, then can identify a specific projection matrix with the Gaussian distribution. Overall, theoretical guarantees and empirical results demonstrate that this method efficiently achieves a deterministic solution that satisfies the JL guarantee.
strength_weakness: strength_weakness: It would be helpful to understand the main results if section 4 could be more organized, such as using subsections.
questions: questions: 1. Could you explain the statement in lines 215-216? Are the values of 1/(3n) and 1/3 derived based on the chosen value of $\epsilon$?

2. Line 336, there is a typo  “appplicability”. 

3. Regarding the notation $\rho$-second-order stationary points($\rho$-SOSP), the paper uses $\rho$-second-order stationary points in some sections and uses $\rho$-SOSP in others.
limitations: limitations: No potential negative societal impact.
paper_summary: paper_summary: The paper considers using the optimization method to "learn" the Johnson Lindenstrauss transform. The paper first shows that the naive objective may not be good enough -- there are stationary points that are sub-optimal. Instead, they consider the way that optimize the random Gaussian space rather than the concrete matrix. Then the authors give an optimization method and show that using this way, every stationary point is a good point and claim that this gives way to deterministic learns the JL matrix. Finally, the paper gives an experiment that shows the advantages of the proposed method.
strength_weakness: strength_weakness: I am still confused about some parts of the paper. I can raise the score if the authors can adjust this. (see the below question)
questions: questions: 1. I am still a little confused about the main conclusion of the results of this paper. That is -- it gives a deterministic construction of the JL lemma, or it gives a better optimization way and it works well empirically? (as the author mentioned, the bound of the JL lemma can not be improved)
 
2. The equation (4) is about probability, and B.1 says they use Monte Carlo sampling, however, would it means that the proposed method still contains some randomness part? 

3. In the experiment section, the paper compares the proposed method with the JL lemma. It will make this stronger if the comparison with equation (1) is also given.
limitations: limitations: See above questions.
paper_summary: paper_summary: The paper proposes to calculate the embedding matrices used in the statement of the Johnson-Lindenstrauss lemma using optimization instead of randomization. The proposed algorithm is a Hessian descent. Authors prove that the algorithm finds the matrix of minimum distortion. Numerical results display the findings
strength_weakness: strength_weakness: The paper lacks complexity analysis of the algorithm. The algorithm proposed requires a full eigenvalue decomposition at every step. It is prohibitive to use this method in any practical scenario. A discussion on the complexity and how to scale the method up (using randomized methods??) would be nice.
questions: questions: The paper's main claim is that mean distortion embeddings are computationally well studied: spectral methods (SVD / eigenvalue) methods calculate those. Authors claim that the min (instead of mean) distortion method is what they want to find. Can authors explain why the relaxation of f* to f using the probability bound in Eq (2) and invoking union bound to obtain (3) does not reduce the max to a sum. The entire promise of the method is to work with the max directly to minimize distortion. It appears that relaxation of max to a sum using the union bound drops the nice property which was the primary motivation of the work. Can you explain what I am missing of misinterpreting?
limitations: limitations: scaling the method up and a clear argument why relaxation of Eq (3) is not making the problem trivial would impove the paper.
paper_summary: paper_summary: This work shows that a deterministic optimization procedure can find a matrix $A$ that satisfies the Johnson Lindenstrauss guarantee. That is, a matrix $A$ maps a set of $n$ vectors to a lower dimensional space while preserving all pairwise distances up to some chosen multiplicative distortion. Typically, $A$ is constructed by sampling it from a random matrix distribution with i.i.d. entries. The authors prove that attempting to directly optimize the entries of $A$ through an optimization procedure by minimizing the maximum distortion is prone to being stuck at local minima. However, the authors show that by optimizing the mean of each entry and the entry-wise variance of the distribution $A$ is sampled from, one can maintain a fixed probability of $A$ being a JL-embedding while at the same time guaranteeing the entry-wise variance ‘sigma’ goes to zero. They then show that, when ‘sigma’ is sufficiently small, one may use the optimized expectation of $A$ as the embedding matrix while only slight increasing the maximum distortion, thereby deterministically finding the desired JL embedding matrix $A$. They show that $\rho$-SOSPs (second order stationary points) have sufficiently low variance when $\rho$ is small, and finally show that a method for finding $\rho$-SOSPs suffices to solve the designed optimization problem.
strength_weakness: strength_weakness: The iterative method to find the matrix $A$ takes $\operatorname{poly}(n, k, d)$ steps, i.e., the complexity is proven to be polynomial but not explicitly determined. Since the paper is primarily theoretical with only limited experiments, it is unclear how efficient this method is in practice.

While the results seem very interesting theoretically, the paper could be strengthened by pointing out some practical applications where this improved deterministic JL embedding would be useful. In the applications I am familiar with, oblivious JL embeddings are needed due to the large number of points in the high-dimensional space (e.g., preserving k-means loss). The authors point to embeddings in deep learning as motivation. It is unclear to me as to how the authors expect progress in understanding deterministic JL embeddings to relate to these embeddings in deep learning. Additional clarification of this point would be helpful.
questions: questions: In the conclusion, you mention the potential for this approach in applications beyond the Johnson Lindenstrauss setting. In your approach for the JL setting, you upper bound the failure probability of the distortion guarantee via the union bound in eqn. (3). This formulation of the objective function seems difficult to translate to other sketching guarantees on $A$ (e.g., projection cost preservation, L2-subspace embedding, affine embedding). Is there any intuitive reason why it may be possible to formulate a relaxed differentiable objective function when the embedding guarantee must hold over an uncountable number of points?

How does learning a JL embedding relate to learning embeddings for application areas discussed in the introduction? In particular, how do you see the results of this paper affecting that line of work? As mentioned above, I think it would be helpful to expand on the link between your result and the motivation of deep learning embeddings given in the intro.
limitations: limitations: Yes

Comments: comment: Thank you for your response.

Regarding the first point, it's important to note that the calculations of the gradient involve **only** cumulative distribution functions of known distributions, i.e. the chi-squared distribution, which simplifies the computations. We acknowledge that this can be explained more clearly and we will incorporate a discussion on the computational aspects in our work.

We appreciate all your valuable feedback and will incorporate it to strengthen our paper. Thank you for your insightful comments and suggestions.
**Q: [...] or is it theoretically proven and related to the objective function that we are optimizing?**

A: It is related to the objective function and theoretically proven, so changing the learning rate won’t resolve it. More generally, to avoid getting stuck at any point during the optimization process, we require either a second-order method or a first-order method with some **randomness**. We expand further on this in the next question.

**Q: Can you share a reference or a more formal statement on this?**

A: Our main result is that the second-order stationary points of Equation (4) satisfy the JL guarantee (Theorem 2). To achieve this, we employ a deterministic second-order method, as detailed in our paper, and thus we achieve a **truly derandomized construction** of a JL matrix.

Additionally, it has been shown that first-order methods with randomness also converge to second-order stationary points with high probability, see for example the work of Jin et al. (2017) on saddle points. Thus, whether using a deterministic second-order method or a first-order method with some degree of randomness, we can reach a second-order stationary point, which is precisely what is required to ensure the JL guarantee.

We hope this clarifies your questions.
can you explain this? 
> second-order optimization is necessary because, without it, using first-order optimization without any randomization results in the optimization getting stuck at the initial point."

is this an empirical finding on the numerical experiments you ran? did you try a larger learning rate? 

or is it theoretically proven and related to the objective function that we are optimizing?



> Generally, to avoid getting stuck throughout the optimization process, you either need second-order information or first-order information combined with randomization. Stochasticity ensures that progress can always be made.

Can you share a reference or a more formal statement on this?
Thanks for the response. We answer your follow-up questions below.

**Q: is there a reason why a first-order optimizer would not work and we actually need a second-order (more expensive) optimizer? First order optimizers are not necessarily stochastic / randomized.**

A simple answer to this is that second-order optimization is necessary because, without it, using first-order optimization **without any randomization** results in the optimization getting stuck at the initial point. Generally, to avoid getting stuck throughout the optimization process, you either need second-order information or first-order information combined with randomization. Stochasticity ensures that progress can always be made. 

**Q: If we remove the $\sigma^2$ term from Equation (4), what estimator do we get? is it not SVD?**

If we remove the $\sigma^2$ term from Equation (4), we would not obtain SVD. Equation (4) provides a bound on the probability of generating a matrix that achieves **worst-case distortion** larger than a specified threshold $\epsilon$. In contrast, PCA, which relies on SVD, finds an embedding that minimizes **average distortion** of distances between the original and the embedded vectors.

By fixing (or as you mentioned, removing) $\sigma^2$, we would essentially be looking to optimize the mean matrix $M$ for that specific $\sigma^2$. However, completely removing $\sigma^2$ from the optimization process would create an undesirable situation. While we might obtain a “better” mean matrix $M$, we would lose the intended derandomization at the end of the optimization, effectively resulting in $\sigma^2 = 0$.
Thanks for the response. I will maintain my score.
Thanks for the response.

- Derandomization: Thanks for pointing this out. However, the calculation of the closed form of expression (4) seems to be non-trivial and some discussion about the computation should be included to support the claim (e.g., the time complexity). This is the main reason I keep my score. 

- Experiment: I agree that "optimizing Equation 1 directly in the matrix space can lead to bad solutions with high distortion" However, it was shown in the specific hard instance and the comparison over the real-world datasets will still be interesting.
thanks for the rebuttal. Follow-up questions: 

* is there a reason why a first-order optimizer would not work and we actually need a second-order (more expensive) optimizer? First order optimizers are not necessarily stochastic / randomized. 

* If we remove the $\sigma^2$ term from Equation (4), what estimator do we get? is it not SVD?
Thanks to the authors for their reply. I will maintain my score.

> Is there any intuitive reason why it may be possible to formulate a relaxed differentiable objective function when the embedding guarantee must hold over an uncountable number of points?

I agree that this paper is valuable for raising the possibility of the general approach you've described being used to achieve deterministic optimization based algorithms for these other sketching guarantees. However, doing so indeed seems quite challenging at first glance.

> How does learning a JL embedding relate to learning embeddings for application areas discussed in the introduction? In particular, how do you see the results of this paper affecting that line of work? 

This is a very interesting interpretation of the results in relation to deep learning. Thanks for clarifying.