Convex Relaxation for Solving Large-Margin Classifiers in Hyperbolic Space

Introduction

The d-dimensional hyperbolic space H d is the unique simply-connected Riemannian manifold with a constant negative sectional curvature -1. Its exponential volume growth with respect to radius motivates representation learning of hierarchical data using the hyperbolic space. Representations embedded in the hyperbolic spaces have demonstrated significant improvements over their Euclidean counterparts across a variety of datasets, including images (Khrulkov et al., 2020), natural languages (Nickel & Kiela, 2017), and complex tabular data such as single-cell sequencing (Klimovskaia et al., 2020).

On the other hand, learning and optimization on hyperbolic spaces are typically more involved than that on Euclidean spaces. Problems that are convex in Euclidean spaces become constrained non-convex problems in hyperbolic spaces. The hyperbolic Support Vector Machine (HSVM), as explored in recent studies Cho et al. (2019); Chien et al. (2021), exemplifies such challenges by presenting as a non-convex constrained programming problem that has been solved predominantly based on projected gradient descent. Attempts have been made to alleviate its non-convex nature through reparametrization (Mishne et al., 2023) or developing a hyperbolic perceptron algorithm that converges to a separator with finetuning using adversarial samples to approximate the large-margin solution (Weber et al., 2020). To our best knowledge, these attempts are 2020), and medical imaging (Wang & Singer, 2013). Some results on the tightness of SDP have been analyzed on a per-problem basis (Bandeira et al., 2017;Brynte et al., 2022;Zhang, 2020). On the other hand, moment-sum-of-squares relaxation, originated from algebraic geometry (Blekherman et al., 2012;Lasserre, 2001), has been studied extensively from a theoretical perspective and has been applied for certifying positivity of functions in a bounded domain (Henrion & Lasserre, 2005). Synthesizing the work done in the control and algebraic geometry literature and geometric machine learning works is under-explored.

Preliminaries

Hyperbolic Space (Lorentz Manifold): define Minkowski product of two vectors x, y → R d+1 as x ↑ y = x 0 y 0 ↓ d i=1 x i y i . A d-dimensional hyperbolic space (Lorentz formulation) is a submanifold embedded in R d+1 defined by,

Tangent Space: a tangent space to a manifold at a given point x → H d is the local linear subspace approximation to the manifold, denoted T x H d . In this case the tangent space is a Euclidean vector space of dimension d written as

Exponential & Logarithmic Map: the exponential map exp x (.) : T x H d ↔ H d is a transformation that sends vectors in the tangent space to the manifold. The logarithmic map log x (.) :

Figure 1: Straight line (red) on Lorentz manifold H 2 as the intersection between a hyperplane and the manifold, presented similarly in Cho et al. (2019).

Exponential and logarithmic maps serve as bridges between Euclidean and hyperbolic spaces, enabling the transfer of notion, such as distances and probability distributions, between these spaces. One way is to consider Euclidean features as residing within the tangent space of the hyperbolic manifold's origin.

From this standpoint, distributions on hyperbolic space can be obtained through exp 0 .

Hyperbolic Decision Boundary: straight lines in the hyperbolic space are intersections between d-dimensional hyperplanes passing through the origin and the manifold H d . Suppose w → R d+1 is the normal direction of the plane, then the plane and hyperbolic manifold intersect if and only if w↑w < 0.

From this viewpoint, each straight line in the hyperbolic space can be parameterized by w and can be considered a linear separator for hyperbolic embeddings. Hence, we can define a decision function h w (.), by the Minkowski product of the feature with the decision plane, as the following,

where w ↑ = [↓w 0 , w 1 , ..., w d ]. A visualization is presented in Figure 1.

Stereographic Projection: we visualize H 2 by projecting Lorentz features isometrically to the Poincaré space B d . Denote Lorentz features as x = [x 0 , x 1 , ..., x d ], then its projection is given by x = [ x1 1+x0 , ...,

Decision boundaries on the Lorentz manifold are mapped to arcs in the Poincaré space. The proof is deferred to Appendix A.2. Cho et al. (2019) proposed the hyperbolic support vector machine which finds a max-margin separator where margin is defined as the hyperbolic point to line distance. We demonstrate our results in a binary classification setting. Extension to multi-class classification is straightforward using Platt-scaling (Platt et al., 1999) in the one-vs-rest scheme or majority voting in one-vs-one setting.

Original Formulation of the HSVM

Suppose we are given {(x i , y i ) :

The hard-margin HSVM is formulated as,

whereas the soft-margin version allows misclassification using

where G → R (d+1)↔(d+1) is a diagonal matrix with diagonal elements diag(G) = [↓1, 1, 1, ..., 1] (i.e. all ones but the first being -1), to represent the Minkowski product in a Euclidean matrix-vector product manner and is the source of indefiniteness of the problem. In the soft-margin case, the hyperparameter C ↭ 0 controls the strength of penalizing misclassification. This penalty scales with hyperbolic distances, defined by l(z) = max(0, arcsinh(1) ↓ arcsinh(z)).

As C approaches infinity, we recover the hard-margin formulation from the soft-margin one. In the rest of the paper we focus on analyzing relaxations to the soft-margin formulation in Equation ( 6) as these relaxations can be applied to both hyperbolic-linearly separable or unseparable data.

To solve the problem e"ciently, we have two observations that lead to two adjustments in our approach. Firstly, although the constraint involving w is initially posited as a strict inequality, practical considerations allow for a relaxation. Specifically, when equality is achieved, w T Gw = 0, the separator is not on the manifold and assigns the same label to all data samples. However, with su"cient samples for each class in the training set and an appropriate regularization constant C, the solver is unlikely to default to such a trivial solution. Therefore, we may substitute the strict inequality with a non-strict one during implementation. Secondly, the penalization function, l, is not a polynomial. Although projected gradient descent is able to tackle non-polynomial terms in the loss function, solvers typically only accommodate constraints and objectives expressed as polynomials. We thus take a Taylor expansion of the arcsinh term to the first order so that every term in the formulation is a polynomial. This also helps with constructing our semidefinite and moment-sum-of-squares relaxations later on, which is presented in Appendix A.1 in detail. The new formulation of the soft-margin HSVM outlined in Equation ( 6) is then given by,

where ω i for i → [n] are the slack variables. More specifically, given a sample (x i , y i ), if ω i = 0, the sample has been classified correctly with a large margin; if ω i → (0, 1 ↗ 2 ], the sample falls into the right region but with a small hyperbolic margin; and if ω i > 1 ↗ 2 , the sample sits in the wrong side of the separator. We defer a detailed derivation of Equation ( 7) to Appendix A.1 and its extension to curvatures other than ↓1 in Appendix A.3.

Semidefinite Formulation

Note that Equation ( 7) is a non-convex quadratically-constrained quadratic programming (QCQP) problem, we can apply a semidefinite relaxation (SDP) (Shor, 1987). The SDP formulation is given by

where decision variables are highlighted in bold and that the last constraint stipulates the concatenated matrix being positive semidefinite, which is equivalent to W ↓ ww T ⇒ 0 by Schur's complement lemma. In this SDP relaxation, all constraints and the objective become linear in (W , w, ω), which could be easily solved. Note that if additionally we mandate W to be rank 1, then this formulation would be equivalent to Equation (7) or otherwise a relaxation. Moreover, it is important to note that this SDP does not directly yield decision boundaries. Instead, we need to extract w ↘ from the solutions (W , w, ω) obtained from Equation (8). A detailed discussion of the extraction methods is deferred to Appendix B.1.

Moment-Sum-of-Squares Relaxation

The SDP relaxation in Equation (8) may not be tight, particularly when the resulting W has a rank much larger than 1. Indeed, we often find W to be full-rank empirically. In such cases, moment-sum-ofsquares relaxation may be beneficial. Specifically, it can certifiably find the global optima, provided that the solution exhibits a special structure, known as the flat-extension property (Curto & Fialkow, 2005;Henrion & Lasserre, 2005).

We begin by introducing some necessary notions, with a more comprehensive introduction available in Appendix C. We define the relaxation order as ε ↭ 1 and our decision variables as q = (w, ω) → R n+d+1 . Our objective, p(q), is a polynomial of degree 2ε with input q, where its coe"cient is defined such that p(q) = 1 2 w T Gw+C n i=1 ω i , thus matching the original objective. Hence, the polynomial p(•) has s(m, 2ε) := m+2ε 2ε number of coe"cients, where m = n + d + 1 is the dimension of decision variables. Additionally, we define z → R s (m,2ε) as the Truncated Multi-Sequence (TMS) of degree 2ε, and we denote a linear functional f associated with this sequence as

which is the inner product between the coe"cients of polynomial p and the vector or real numbers z. The vector of monomials up to degree ε generated by q is denoted as [q] ε . With all these notions established, we can then define the moment matrix of ε-th degree, M ε [z], and localizing matrix of ε-th degree for polynomial g, L ε,g [z], as the followings,

where s is the max degree such that 2s + deg(g) ↫ 2ε, [q] ε [q] T ε is a matrix of polynomials with size s(m, ε) by s(m, ε), and all the inner products are applied element-wise above. For example, if n = 1 and d = 2 (i.e. 1 data sample from a 2-dimensional hyperbolic space), the degree-2 monomials generated by q = (w 0 , w 1 , w 2 , ω 1 ) are

With all these definitions established, we can present the moment-sum-of-squares relaxation (Nie, 2023) to the HSVM problem, outlined in Equation ( 7), as (Moment) min z↓R s(m,2ω) ⇑vec(p), z⇓ .

Note that g(q) ↭ 0, as previously defined, serves as constraints in the original formulation. Additionally, when forming the moment matrix, the degree of generated monomials is s = ε ↓ 1, since all constraints in Equation ( 7) has maximum degree 1. Consequently, Equation ( 13) is a convex programming and can be implemented as a standard SDP problem using mainstream solvers. We further emphasize that by progressively increasing the relaxation order ε, we can find increasingly better solutions theoretically, as suggested by Lasserre (2018). However, moment-sum-of-squares relaxation does not scale with the data size due to the combinatorial factors in the dimension of truncated multi-sequence z, leading to prohibitively slow runtimes and excessive memory consumption. To address this issue, we exploit the sparsity pattern inherent in this problem: many generated monomial terms do not appear in the objective or constraints. For instance, there is no cross-terms among the slack variables, such as ω i ω j for i ⇔ = j → [n]. Specifically, in this problem, we observe a star-shaped sparsity structure, as ilustrated in Figure 2. We observe that, by defining sparsity groups as q (i) = (w, ω i ), two nice structural properties can be found: first, the objective function involves all the sparsity groups, {q (i) } n i=1 and can be decomposed as sum of polynomials p (i) , which involves only q (i) and thus has a smaller number of coe"cients s(m ↑ , 2ε), where m ↑ = d + 2. Second, each constraint is exclusively associated with a single group q (i) for a specific i. For the remaining constraint, w T Gw, we could assign it to group i = 1 without loss of generality. Hence, by leveraging this sparsity property, we can reformulate the moment-sum-of-squares relaxation equivalently into its sparse version,

where B is an index set of the moment matrix to entries generated by w along, ensuring that each moment matrix with overlapping regions share the same values as required. We refer the last constraint as the sparsebinding constraint. The nice thing about Equation ( 14) is that it is equivalent to Equation ( 13) and has a fewer number of decision variables.

Unfortunately, our solution empirically does not satisfy the flat-extension property and we cannot not certify global optimality. Nonetheless, in practice, it achieves significant performance improvements in selected datasets over both projected gradient descent and the SDP-relaxed formulation. Similarly, this formulation does not directly yield decision boundaries and we defer discussions on the extraction methods to Appendix B.2.

Experiments

We validate the performances of semidefinite relaxation (SDP) and sparse moment-sum-of-squares relaxations (Moment) by comparing various metrics with that of projected gradient descent (PGD) on a combination of synthetic and real datasets. The PGD implementation follows from adapting the MATLAB code in Cho et al. (2019), with learning rate 0.001 and 2000 epochs for synthetic and 4000 epochs for real dataset and warm-started with a Euclidean SVM solution. Datasets. For synthetic datasets, we construct Gaussian and tree embedding datasets following Cho et al. (2019); Mishne et al. (2023); Weber et al. (2020). Regarding real datasets, our experiments include two machine learning benchmark datasets, CIFAR-10 Krizhevsky et al. (2009) and Fashion-MNIST Xiao et al. (2017) with their hyperbolic embeddings obtained through standard hyperbolic embedding procedure (Chien et al., 2021;Khrulkov et al., 2020;Klimovskaia et al., 2020) to assess image classification performance. Additionally, we incorporate three graph embedding datasets, such as football, karate, and polbooks obtained from Chien et al. (2021), to evaluate the e!ectiveness of our methods on graph-structured data. We also explore cell embedding datasets, including Paul Myeloid Progenitors developmental dataset (Paul et al., 2015), Olsson Single-Cell RNA sequencing dataset (Olsson et al., 2016), Krumsiek Simulated Myeloid Progenitors dataset (Krumsiek et al., 2011), and Moignard blood cell developmental trace dataset from single-cell gene expression (Moignard et al., 2015), where the inherent geometry structures well fit into our methods.

We emphasize that all features are on the Lorentz manifold, but visualized in Poincaré manifold through stereographic projection if the dimension is 2.

Evaluation Metrics. The primary metrics for assessing model performance are average training and testing loss, accuracy, and weighted F1 score under a stratified 5-fold train-test split scheme. Furthermore, to assess the tightness of the relaxations, we examine the relative suboptimality gap, defined as

where f ↘ is the unknown optimal objective value, p ↘ is the objective value of the relaxed formulation, and f is the objective associated to the max-margin solution recovered from the relaxed model. Clearly p ↘ ↫ f ↘ ↫ f , so if ϑ ↖ 0, we can certify the exactness of the relaxed model.

Implementations Details. We use MOSEK (ApS, 2022) in Python as our optimization solver without any intermediate parser, since directly interacting with solvers save substantial runtime in parsing the problem.

MOSEK uses interior point method to update parameters inside the feasible region without projections. All experiments are run and timed on a machine with 8 Intel Broadwell/Ice Lake CPUs and 40GB of memory.

Results over multiple random seeds have been gathered and reported.

We first present the results on synthetic Gaussian and tree embedding datasets in Section 4.1, followed by results on various real datasets in Section 4.2.

Synthetic Dataset

Synthetic Gaussian. To generate a Gaussian dataset on H d , we first generate Euclidean features in R d and lift to hyperbolic space through exponential map at the origin, exp 0 , as outlined in Equation (3). We adjust the number of classes K → {2, 3, 5} and the variance of the isotropic Gaussian s → {0.4, 0.6, 0.8, 1.0}.

Published in Transactions on Machine Learning Research (04/2025) Three Gaussian embeddings in d = 2 are selected and visualized in Figure 3 and performances with C = 10 for the three dataset are summarized in Table 1.

In general, we observe a small gain in average test accuracy and weighted F1 score from SDP and Moment relative to PGD. Notably, we observe that Moment often shows more consistent improvements compared to SDP, across most of the configurations. In addition, Moment gives smaller optimality gaps ϑ than SDP. This matches our expectation that Moment is tighter than the SDP.

Although in some case, for example when K = 5, Moment achieves significantly smaller losses compared to both PGD and SDP, it is generally not the case. We emphasize that these losses are not direct measurements of the max-margin hyperbolic separators' generalizability; rather, they are combinations of margin maximization and penalization for misclassification that scales with C. Hence, the observation that the performance in test accuracy and weighted F1 score is better, even though the loss computed using extracted solutions from SDP and Moment is sometimes higher than that from PGD, might be due to the complicated loss landscape. More specifically, the observed increases in loss can be attributed to the intricacies of the landscape rather than the e!ectiveness of the optimization methods. Based on the accuracy and F1 score results, empirically SDP and Moment methods identify solutions that generalize better than those obtained by running gradient descent alone. We provide a more detailed analysis on the e!ect of hyperparameters in Appendix E.2 and runtime in Table 5. Decision boundary for Gaussian 1 is visualized in Figure 5 as an example. Synthetic Tree Embedding. As hyperbolic spaces are good for embedding trees, we generate random tree graphs and embed them to H 2 following Mishne et al. ( 2023). Specifically, we label nodes as positive if they are children of a specified node and negative otherwise. Our models are then evaluated for subtree classification, aiming to identify a boundary that includes all the children nodes within the same subtree.

Such task has various practical applications. For example, if the tree represents a set of tokens, the decision boundary can highlight semantic regions in the hyperbolic space that correspond to the subtrees of the data graph. We emphasize that a common feature in such subtree classification task is data imbalance, which usually lead to poor generalizability. Hence, we aim to use this task to assess our methods' performances under this challenging setting. Three embeddings are selected and visualized in Figure 3 and performance is summarized in Table 1. The runtime of the selected trees can be found in Table 5. Decision boundary of tree 2 is visualized in Figure 6.

Similar to the results of synthetic Gaussian datsets, we observe better performance from SDP and Moment compared to PGD, and due to data imbalance that gradient-based methods typically struggle with, we have a larger gain in weighted F1 score in this case. In addition, we observe large optimality gaps for SDP but very tight gap for Moment, certifying the optimality of Moment even when class-imbalance is severe.

Real Dataset

Real datasets consist of embedding of various sizes and number of classes in H 2 , visualized in Figure 4. We first report performances of three models using one-vs-rest training scheme, described in Appendix D, in Tables 6 to 8 for C → {0.1, 1.0, 10} respectively, and report aggregated performances, by selecting the one with the highest average test weighted F1 score, in Table 2. In general, we observe that Moment achieves the best test accuracy and weighted F1 score, particularly in biological datasets with clear hyperbolic structures, and have smaller optimality gaps compared to SDP relaxation, for nearly all selected data. However, it is important to note that the optimality gaps of these two methods remain distant from zero, suggesting that these relaxations are not tight enough for these datasets. Nevertheless, both relaxed models significantly outperform projected gradient descent (PGD) by a wide margin. Furthermore, our observations reveal that in the one-vs-rest training scheme, PGD shows considerable sensitivity to the choice of the regularization parameter C from Tables 6 to 8, whereas SDP and Moment are less a!ected, demonstrating better stability and consistency across di!erent C's.

One critical drawback of semidefinite and sparse moment-sum-of-squares relaxation is that they do not scale e"ciently with an increase in data samples, resulting in excessive consumption of time and memory, for example, CIFAR10 and Fashion-MNIST using a one-vs-rest training scheme. The workaround is one-vsone training scheme, where we train for O(K 2 ) number of classifiers among data from each pair of classes and make final prediction decision using majority voting. We summarize the performance in Table 3 by aggregating results for di!erent C in Tables 9 to 11 as in the one-vs-rest case. We observe that in one-vsone training, the improvement in general from the relaxation is not as significant as it in the one-vs-rest scheme, and SDP relaxation now gives the best performance in average test accuracy and test F1, albeit with large optimality gaps. Note that in the one-vs-one scheme, PGD is more consistent across di!erent C's, potentially because each subproblem-binary classifying one class against another-contains less data compared to one-vs-rest, making it easier to identify solutions.

A more detailed analysis on the e!ect of regularization C and runtime comparisons are provided in Appendix E.3 and Table 12. Table 2: Real dataset 5-fold test accuracy, F1, and optimality gap with one-vs-rest training. Best metrics based on weighted F1 is reported here aggregating from Tables 6 to 8. data test acc test f1 (micro) ϑ PGD SDP Moment PGD SDP Moment SDP Moment football 40.87% ± 4.43% 32.17% ± 4.43% 37.39% ± 4.43% 0.29 ± 0.03 0.23 ± 0.04 0.26 ± 0.03 0.3430 0.0999 karate 50.00% ± 6.39% 50.00% ± 6.39% 50.00% ± 6.39% 0.66 ± 0.06 0.66 ± 0.06 0.66 ± 0.06 0.9155 0.0818 polbooks 84.76% ± 1.90% 84.76% ± 1.90% 84.76% ± 1.90% 0.79 ± 0.02 0.80 ± 0.03 0.80 ± 0.03 0.1711 0.0991 krumsiek 81.78% ± 2.66% 82.56% ± 2.01% 86.47% ± 0.64% 0.79 ± 0.03 0.80 ± 0.03 0.84 ± 0.00 0.7519 0.0921 moignard 63.37% ± 0.70% 63.68% ± 1.75% 63.78% ± 1.57% 0.62 ± 0.01 0.60 ± 0.02 0.60 ± 0.02 0.0325 0.0396 olsson 74.27% ± 5.15% 79.63% ± 3.54% 81.20% ± 3.68% 0.69 ± 0.07 0.77 ± 0.04 0.79 ± 0.04 0.4118 0.0976 paul 54.85% ± 1.26% 53.72% ± 2.42% 64.71% ± 2.36% 0.48 ± 0.02 0.47 ± 0.03 0.61 ± 0.02 0.4477 0.0861 myeloidprogenitors 69.34% ± 3.81% 70.12% ± 3.28% 76.84% ± 2.04% 0.66 ± 0.05 0.67 ± 0.04 0.75 ± 0.02 0.6503 0.1074

Limitations

Although our proposed SDP and sparse moment relaxation methods achieve generally better classification performance compared to the PGD baseline in the one-vs-rest setting -and slightly less so in the one-vs-one setting -the computational cost is significant due to their longer runtimes, as previously mentioned. One naturally suspects that given equivalent computational resources (including time and memory), PGD might achieve performance comparable to the sparse moment relaxation methods. This scenario becomes especially plausible considering the large hyperparameter space available for PGD: a more aggressive grid search on learning rates, total number of epochs, regularization parameters, and initializations could potentially yield results that rival or even surpass those of the moment-based method.

To explore this further, we provide an additional comparison between our methods and PGD under approximately equal computational time constraints using the one-vs-rest scheme. Specifically, according to Table 5, the average runtime per training fold for the Moment method is approximately two orders of magnitude larger than that for PGD. Therefore, we randomly generated 100 di!erent initializations for PGD, with each initialization vector drawn uniformly from the unit sphere subject to the condition w T Gw > 0, and report the best performance among these random initializations in Table 4 below. It's worth noting that for PGD, we follow Cho et al. (2019), where the initialization w is determined based on the solution w ↑ of a soft-margin SVM in the ambient Euclidean space of the Lorentz model. This choice typically ensures better initialization, contributing to the stability of optimization, especially when numerous local optima may exist. Consequently, we observed that the best PGD performance among 100 di!erent random initializations do not significantly improve its performance in this regime. Nevertheless, we do not exclude the possibility that more strategic hyperparameter tuning under similar computational constraints might yield improved PGD performance

Based on these insights, we o!er practical recommendations for method selection tailored to dataset characteristics and scale. For small-scale, low-dimensional hyperbolic datasets, particularly in 2D, where a hyperparameter-free solution is desired, the SDP and sparse moment relaxation methods clearly excel. Additionally, we advocate exhibiting significant class overlap. For instance, arguably Gaussian 1 and 3 have greater overlap compared to Gaussian 2 (Figure 3), leading to a notable performance drop in the PGD method. We suspect that more mixing likely leads PGD to poor local minima, although this conjecture warrants further quantitative investigation. For larger dataset such as MNIST or comparable modern machine learning benchmarks, in which generalization to unseen data takes priority over strict global optimality, we recommend PGD due to its significantly lower computational demands compared to our relaxed proposals.

Discussions

In this paper, we provide a stronger classification performance on hyperbolic support vector machine using semidefinite and sparse moment-sum-of-squares relaxations on the hyperbolic support vector machine problem compared to projected gradient descent. In the one-vs-rest settings, we observe that they achieve better classification accuracy and F1 score than the existing PGD approach on both simulated and real dataset. Additionally, we discover small optimality gaps for moment-sum-of-squares relaxation, which approximately certifies global optimality of the moment solutions.

Perhaps the most critical drawback of SDP and sparse moment-sum-of-squares relaxations is their limited scalability. The runtime and memory consumption grows quickly with data size and dimension. We attempted dividing this challenge into sub-tasks, such as using one-vs-one training scheme, to cut overall runtime and memory, but observe a smaller performance advantage over the PGD baseline as they seem to perform better in this alternative setting. Also, if given the same computational resources as moment-sum-Published in Transactions on Machine Learning Research (04/2025) of-squares, we suspect that PGD may be able to perform on par with the moment-relaxation method using a rather greedy hyperparameter hunt such as specifically tunning the learning rates and regularization factors.

But it is not completely without hope. For relatively large datasets, we could develop more heuristic approaches for solving our relaxed optimization problems to achieve runtime comparable with projected gradient descent. Combining the GD dynamic with interior point iterates in a problem-dependent manner could be useful Yang et al. (2023).

It remains to show if we have performance gain in either runtime or optimality by going through the dual of the problem or by designing feature kernels that map hyperbolic features to another set of hyperbolic features (Lensink et al., 2022). Nonetheless, we believe that our work introduces a valuable perspectiveapplying SDP and Moment relaxations -to the geometric machine learning community.

By introducing the slack variable ω i as the penalty term in Equation ( 6), we can rewrite Equation ( 6) into

then by rearranging terms and taking sinh on both sides, it follows that

where the last equality follows from hyperbolic trig identities. To make it ready for moment-sum-of-squares relaxation, we turn the function into a polynomial constraint by taking the taylor expansion up to some odd orders (we need monotonic decreasing approximation to g(.), so we need odd orders).

If taking up to the first order, we relax the problem into Equation ( 7)foot_0 . If taking up to the third order, we relax the original problem to,

It's worth mentioning that we expect the lower bound gets tighter as we increase the order of Taylor expansion. However, once we apply the third order Taylor expansion, the constraint is no longer quadratic, eliminating the possibility of deriving a semidefinite relaxation. Instead, we must rely on moment-sum-ofsquares relaxation, potentially requiring a higher order of relaxation, which may be highly time-costly.

It is also worth noting that such a Taylor expansion is consistent with the arcsinh when we allow the curvature of the hyperbolic space to vary. Specifically, one can show that as the curvature approaches 0 (i.e. as the hyperbolic space tends to Euclidean), both the arcsinh formulation in Equation ( 6) and Equation ( 7) formulation recovers the standard Euclidean SVM problem. This further justifies the naturalness of taking such a Taylor expansion.

A.3 An extension to curvature other than -1

Throughout the paper, we assume that the hyperbolic space has a constant negative curvature ↓1 without loss of generality. Our analysis and relaxation easily extend to curvature other than ↓1.

Define c as the negative curvature. That is, all the discussions in the main text are about c = 1. To derive the prooblem formulation and associated relaxation, we can follow Table 1 in Skopek et al. (2019) and retrace the proof in the supplemental material in Chien et al. ( 2021), then we see that the decision function, the objective, and many constraints stay the same except the following slack variable constraint in Equation ( 16):

and with some rearrangements, we have

whose Taylor expansion around ω i = 0 gives

Essentially, this is nothing but adding a curvature factor c into the problem formulation, so that solving the SVM problem in other curvature has no additional theoretical challenges. The rest follows likewise as in Appendix A.1.

The intuition behind why only this single constraint is a!ected can be made as follows: SVM cares about direction of the decision boundary, and if the optimal direction is found for an embedding in one c, the "same" direction is also optimal for the same embedding stretched to other C; but since the slack variables are associated to distance measures, they need to be scaled accordingly.

B.1 Semidefinite Relaxation

For SDP, we initially construct a set of candidates w derived from (W , w, ω). Then, among candidates in this set, we choose the one that minimizes the loss function in Equation ( 7).

The candidates, denoted as w's, include 1. Scaled top eigendirection: w = ↘ ϖ max u max , where ϖ max and u max are the largest eigenvalue and the eigenvector associated with the largest eigenvaue; 2. Gaussian randomizations: sample w ↙ N (w, W ↓ ww T )foot_2 . We empirically generate 10 samples from this distributions;

3. Scaled matrix columns: if it were the case that W = ww T , then each column of W contains w scaled by some entry within itself. Using columns of W divided by the corresponding entry of w (e.g. divide first column by w 0 , second column by w 1 , and so on), we get d + 1 many candidates w = w's;

4. Nominal solution: w = w, i.e. include w itself as a candidate.

Typically the top eigendirection is selected as the best candidate.

B.2 Moment-Sum-of-Squares Relaxation

In moment-sum-of-squares relaxation, the decision variable is the truncated multi-sequence z, but we could decode the solution from the moment matrix M ε [z] it generates. We are able to extract the part in TMS that corresponds to w = [w 0 , w 1 , ..., w d ], by reading o! these entries from the moment matrix, which is already a good enough solution.

For example, in d = 2, ε = 2, one of the sparcity group, say q (1) consists of [w, ω 1 ], which has monomials generated in Equation ( 12). Define ∝ as a binary operator between two vectors of monomials that generates another vector with monomials given by the unique combinations of the product into vectors, such that

Then, monomials generated can be more succinctly expressed as

and the moment matrix can be expressed in block form as

Note that the value for w (the red part) is contained close to the top left corner of the moment matrix, which provides us good linear separator w in this problem.

C.1 Polynomial Optimization and Dual Cones

Polynomial optimization problem (POP) in the most generic form can be presented as

where p(x) is our polynomial objective and h i (x), g i (x) are our polynomial equality and inequality constraints respectively. However, in general, solving such POP to global optimality is NP-hard (Lasserre, 2001;Nie, 2023). To address this challenge, we leverage methods from algebraic geometry (Blekherman et al., 2012;Nie, 2023), allowing us to approximate global solutions using convex optimization methods.

To start with, we define sum-of-squares (SOS) polynomials as polynomials that could be expressed as a sum of squares of some other polynomials, and we define ![x] to be the collection of SOS polynomials.

More formally, we have

where R[x] denotes the polynomial ring over R.

Next, we recall the definitions of quadratic module and its dual. Given a set of polynomials g = [g 1 , g 2 , ..., g l ], the quadratic module generated by g is defined as

and its degree 2d-truncation is defined as,

where g 0 = 1. It has been shown that the dual cone of Qmodule[g] 2d is exactly the convex cone defined by the PSD conditions of the localizing matrices,

refers to the d th order localizing matrix of g i generated by z, and f z (g i ) = ⇑f z , g i ⇓ = ⇑z, vec(g i )⇓ refers to the linear functional associated with z applied on g i → R[x] 2d . It is worth mentioning that the application of the linear functional f z to the symmetric polynomial matrix g(x)

s is element-wise. Formally speaking, for all g ↑ → Qmodule[g] 2d and for all g → M[g] 2d , we have ⇑g ↑ , g⇓ ′ 0. Similarly, given a set of polynomials h = [h 1 , h 2 , ..., h m ], the ideal generated by h is defined as,

and its degree 2d-truncation is defined as,

where ϖ i 's are also called polynomial multipliers. Interestingly, it is shown that we can perfectly characterize the dual of the sum of ideal and quadratic module,

} is the convex cone defined by the PSD conditions of the localizing matrices.

With these notions setup, we can reformulate the POP above into the following SOS program for arbitrary ε → N as the relaxataion order,

whose optimal value produces a lower bound to p ↘ , i.e. ς ↘ ε △ p ↘ , and its dual problem of the SOS program above is, n,2d) ⇑f z , p⇓ .

This pair of SOS programs is called the moment-sum-of-squares hierarchy first proposed in Lasserre (2001). It is particularly useful as it has been shown that

and {ς ↘ ε } ε and {φε ↘ } ε are two monotonically increasing sequences. In our work, we implement our SOS programs following the dual route.

C.2 Sparse Polynomial Optimization

In this section, we briefly discuss how sparse moment-sum-of-squares is formulated. Using the same sparsity pattern defined in Section 3 (i.e. q (i) = (w, ω i )), we first introduce the notion of correlated sparsity. Definition 1. Correlated Sparsity for an objective p → R[x] and associated set of constraints means 1. For any constraint g i (q), ⇐i → [l], it only involves term in one sparsity group q (i) for some i → [n]

2. The objective can be split into

3. The grouping satisfies the running intesection property (RIP), i.e. for all i → {1, 2, ..., n ↓ 1}, we have

≃ q (s) , for some s ↫ i

In our case, the first property is straightforward. For the second, we may define explicitly p i (q (i) ) = p i (w, ω i ) = 1 2n w T Gw + Cω i so that we get back the original objective after summation. The last property direct follows from the star-shaped structure, i.e. ⇐i → {1, 2, ..., n ↓ 1}, we indeed have ̸ i k=1 q (i) ▽ q (i+1) = w ≃ q (1) . Hence, our sparsity group indeed satisfies all three property and thus we have correlated sparsity in the problem.

With correlated sparsity and data regularity (Putinar's Positivestellentz outlined in Nie ( 2023)), we are able to decompose the Qmodule generated by the entire set of decision variables into the Minkowski sum of Qmodules generated by each sparsity group of variables, e!ectively reducing the number of decision variables in the implementations. For a problem with only inequality constraints, which is our case for HSVM, the sparse POP for our problem reads as max ς ,

and we could derive its dual accordingly and present the SDP form for implementation in Equation ( 14).

D Platt Scaling (Platt et al., 1999)

Platt scaling (Platt et al., 1999) is a common way to calibrate binary predictions to probabilistic predictions in order to generalize binary classification to multiclass classification, which has been widely used along with SVM. The key idea is that once a separator has been trained, an additional logistic regression is fitted on scores of the predictions, which can be interpreted as the closeness to the decision boundary.

In the context of HSVM, suppose w ↘ is the linear separator identified by the solver, then we find two scalars, A, B → R, with

where ↑ refers to the Minkowski product defined in Equation ( 1). The value of A and B are trained on the trained set using logistic regression with some additional empirical smoothing. For one-vs-rest training, we will then have K sets of (A, B) to train, and at the end we classify a sample to the class with the highest probability. See detailed implementation here https://home.work.caltech.edu/ htlin/program/libsvm/doc/platt.py in LIBSVM.

E.1 Visualizing Decision Boundaries

Here we visualize the decision boundary of for PGD, SDP relaxation and sparse moment-sum-of-squares relaxation (Moment) on one fold of the training to provide qualitative judgements.

We first visualize training on the first fold for Gaussian 1 dataset from Figure 3 in Figure 5. We mark the train set with circles and test set with triangles, and color the decision boundary obtained by three methods with di!erent colors. In this case, note that SDP and Moment overlap and give identical decision boundary up to machine precision, but they are di!erent from the decision boundary of PGD method. This slight visual di!erence causes the performance di!erence displayed in Table 1.

We next visualize the decision boundary for tree 2 from Figure 3 in Figure 6. Here the di!erence is dramatic: we visualize both the entire data in the left panel and the zoomed-in one on the right. We indeed observe that the decision boundary from moment-sum-of-squares relaxation have roughly equal distance from points to the grey class and to the green class, while SDP relaxation is suboptimal in that regard but still enclosing the entire grey region. PGD, however, converges to a very poor local minimum that has a very small radius enclosing no data and thus would simply classify all data sample to the same class, since all data falls to one side of the decision boundary. As commented in Section 4, data imbalance is to blame, in which case the final converged solution is very sensitive to the choice of initialization and other hyperparameters such as learning rate. This is in stark contrast with solving problems using the interior point method, where after implementing into MOSEK, we are essentially care-free. From this example, we see that empirically sparse moment-sum-of-squares relaxation finds linear separator of the best quality, particularly in cases where PGD is expected to fail.

E.2 Synthetic Gaussian

To generate mixture of Gaussian in hyperbolic space, we first generate them in Euclidean space, with the center coordinates independently drawn from a standard normal distribution. K such centers are drawn for defining K di!erent classes. Then we sample isotropic Gaussian at respective center with scale s. Finally, we lift the generated Gaussian mixtures to hyperbolic spaces using exp 0 . For simplicity, we only present results for the extreme values: K → {2, 5}, s → {0.4, 1}, and C → {0.1, 10}.

For each method (PGD, SDP, Moment), we compute the train/test accuracy, weighted F1 score, and loss on each of the 5 folds of data for a specific (K, s, C) configuration. We then average these metrics across the 5 folds, for all methods and configurations. To illustrate the performance, we plot the improvements of the average metrics of the Moment and SDP methods compared to PGD as bar plots for 15 di!erent seeds.

Outliers beyond the interquartile range (Q1 and Q3) are excluded for clarity, and a zero horizontal line is marked for reference. Additionally, to compare the Moment and SDP methods, we compute the average optimality gaps similarly, defined in Equation ( 15), and present them as bar plots. Our analysis begins by examining the train/test accuracy and weighted F1 score of the PGD, SDP, and Moment methods across various synthetic Gaussian configurations, as shown in Figures 7 to 10.

Across various configurations, we observe that both the Moment and SDP methods generally show improvements over PGD in terms of train and test accuracy as well as weighted F1 score. Notably, we observe that Moment method often shows more consistent improvements compared to SDP. This consistency is evident across di!erent values of (K, s, C), suggesting that the Moment method is more robust and provide more generalizable decision boundaries. Moreover, we observe that 1. for larger number of classes (i.e. larger K), the Moment method consistently and significantly outperforms both SDP and PGD, highlighting its capability to manage complex class structures e"ciently; and 2. for simpler datasets (with smaller scale s), both Moment and SDP methods generally outperform PGD, where the Moment method particularly shows a promising performance advantage over both PGD and SDP.

Next, we move to examine the train/test loss improvements compared to PGD and optimality gaps comparison across various configurations, shown in Figures 11 to 14. We observe that for K = 5, the Moment method achieves significantly smaller losses compared to both PGD and SDP, which aligns with our previous observations on accuracy and weighted F1 scores. However, for K = 2, the losses of the Moment and SDP methods are generally larger than PGD's. Nevertheless, it is important to note that these losses are not direct measurements of our optimization methods' quality; rather, they measure the quality of the extracted solutions. Therefore, a larger loss does not necessarily imply that our optimization methods are inferior to PGD, as the heuristic extraction methods might significantly impact the loss. Additionally, we observe that the optimality gaps of the Moment method are significantly smaller than those of the SDP method, suggesting that Moment provides better solutions. Interestingly, the optimality gaps of the Moment method also exhibit smaller variance compared to SDP, as indicated by the smaller boxes in the box plots, further supporting the consistency and robustness of the Moment method. Lastly, we compare the computational e"ciency of methods, where we compute the average runtime to finish 1 fold of training for each model on synthetic dataset, shown in Table 5. We observe that sparse moment relaxation typically requires at least one order of magnitude runtime compared to other which to some extent limits the applicability of this method to large scale dataset.

E.3 Real Data

In this section we provide detailed performance breakdown by the choice of regularization C for both onevs-one and one-vs-rest scheme in Tables 6 to 11. In one-vs-rest scheme, we observe that the Moment method consistently outperforms both PGD and SDP across almost all datasets and C in terms of accuracy and F1 scores. Notably, the optimality gaps, ϑ, for Moment are consistently lower than those for SDP, indicating that the Moment method's solution obatin a better gap, which underscore the e!ectiveness of the Moment method in datasets. In one-vs-one scheme however, we observe that the SDP and Moment have comparative performances, both better than PGD. Nevertheless, the optimality gaps of SDP are still significantly larger than the Moment's, for almost all cases.

Similarly, we compare the average runtime to finish 1 fold of training for each model on these real datasets, shown in Table 12. We observe a similar trend: the sparse moment relaxation typically requires at least an order of magnitude more runtime compared to the other methods.

For the implementation in MOSEK, we linearize the l 1 norm term by introducing extra auxiliary variables, which we do not show here. The moment relaxation can be implemented likewise, since this is constraint-wise uncertainty and we preserve the same sparsity pattern so that the same sparse moment relaxation applies.

Remark 1. Since we take a Taylor approximation to first order, Ũx also does not guarantee that all of its elements live strictly on the hyperbolic manifold H d . One could seek the second order approximation to the defining function, which make the robust criterion a quadratic form on the uncertainty parameter z. Remark 2. Equation (33) arises from a l ⇐ uncertainty set. We could of course generalize this analysis to other types of uncertainty sets such as l 2 uncertainty and ellipsoidal uncertainty, each with di!erent number of auxiliary variables introduced in their linearization.