A skewness-aware matrix factorization approach for mesh-structured cloud services

1 A Skewness-aware Matrix Factorization Approach for Mesh-structured Cloud Services Yongquan Fu, Dongsheng Li, Pere Barlet-Ros, Chun Huang, Zhen Huang, Siqi Shen, Huayou Su Online cloud services need to fulfill clients’ requests scalably and fast, otherwise, users’ experience and providers’ revenue could be severely degraded. State-of-the-art cloud services are increasingly deployed as a distributed service mesh. Service to service communication is frequent in the mesh. Unfortunately, problematic events may occur between any pair of nodes in the mesh, therefore, it is vital to maximize the network visibility for efficient network troubleshooting and application optimization. A state-of-the-art approach is to model pairwise RTTs based on a latent factor model represented as a low-rank matrix factorization. A latent factor corresponds to a rank-1 component in the factorization model, and is shared by all node pairs. However, different node pairs usually experience a skewed set of hidden factors like divergent routing paths, dynamic path conditions, or transient network congestions, which cannot be captured by existing latent factor models. In this paper, we propose a skewness-aware matrix factorization method named SMF. We decompose the matrix factorization into the basic units of rank-one latent factors, and progressively combine rank-one factors for different node pairs. We present a unifying framework to automatically and adaptively select the rank-one factors for each node pair, which not only preserves the low rankness of the matrix model, but also adapts to skewed network latency distributions. Over real-world RTT data sets, SMF significantly improves the relative error by a factor of 0.2x to 10x, converges fast and stably, and compactly captures fine-grained local and global network latency structures. I. INTRODUCTION Large-scale cloud services are typically organized as a mesh of micro-services that are deployed over hundreds to thousands of nodes, as illustrated in Figure 1. Service-to-service communication is frequent on the mesh structured service topology. For example, a Web request may traverse thousands of servers to search and aggregate results. A request’s service level agreement (SLA) is determined based on the response from the slowest server. As problematic locations are essentially unpredictable a priori, we need to track RTTs by all nodes and for all nodes [25]. Timely response is vital for ensuring the Quality of Experience (QoE) [49], otherwise, the increased delay significantly affects users’ experience and providers’ revenue [13], [7]. Unfortunately, high latency issues may arise between any node pairs of the service mesh, due to changing routing paths, degraded path conditions, or transient network congestions. Yongquan Fu, Dongsheng Li, Zhen Huang, Siqi Shen, Huayou Su are with the Science and Technology Laboratory of Parallel and Distributed Processing; College of Computer, National University of Defense Technology, Changsha, Hunan province, China. Chun Huang is with College of Computer, National University of Defense Technology, Changsha, Hunan province, China. Pere Barlet-Ros is with UPC BarcelonaTech. Fig. 1. An illustration of a service-mesh scheme [3], [4]. A cloud service consists of a collection of loosely-coupled micro-services hosted on a set of distributed servers. The service-to-service communication is managed by a dedicated service-mesh proxy layer for autonomous traffic control and optimization. The monitoring functionality is located at the proxy layer to collect pairwise RTT status for the service mesh layer. As a result, in order to correctly correlate the network issue with the problematic service, it is vital to tell if a service is affected by a network issue by monitoring pairwise RTTs of the service mesh. Meanwhile, optimizing the performance of the service mesh needs pairwise RTTs. For example, we may choose the detour routing scheme [44], [6], [35] that selects a relay on the direct routing path between a set of hosts and forwards packets via the relay. For example, Akamai SureRoute [40] and Amazon CloudFront content distribution network (CDN) [5] relay packets for end hosts. A straightforward approach is to directly collect all-pair network latency, which requires a quadratic number of probing packets with respect to the system size, which does not scale well. Therefore, in order to provide enough network visibility for troubleshooting, it is of paramount importance to predict missing measurement results. Besides the scaling limitations, the measurement may also disturb normal application traffic, as node resources are usually shared among multiple tenants. To reduce the probing cost, researchers have proposed network-latency prediction methods that embed the service mesh into a low-dimensional vector space and estimate the pairwise RTTs based on the vector distance of the corresponding node pair. The vector-space representation succinctly captures the pairwise network latency matrix, and needs only O(N)probes to predict the full pairwise matrix for Nhosts in edge data centers where PingMesh becomes ineffective. Further, the vector representation can act as the input for various applications [45], [35], [28]. The vector space is usually represented as a low-rank matrix factorization model [37], [33], [34], which represents each node as a vector of variables. These variables are “latent” in the sense that they are not directly observed, but should be inferred via a mathematical model from observed data. A 2 matrix factorization model is well known to be equivalently represented as the sum of a set of rank-one matrices [24]. Given a rank-one matrix represented as the product of two vectors FkGT kfor k≤r, where a constant rdenotes the number of rank-one matrices, this rank-one matrix serves as a “latent factor”, and the vector Fkand Gkserve as the soft memberships of each row and each column of the RTT matrix to this latent factor [15], [47]. As the same set of rankone matrices are shared by all nodes, the matrix factorization model implicitly assumes that all node pairs are affected by the same set of latent factors. However, our analyses in Section III-C show that a real-world network latency matrix is heterogeneous and highly skewed, different node pairs are typically correlated with heterogeneous factors. Consequently, the matrix factorization should be revisited to account for skewed factors. In this paper, we focus on improving the skewness awareness for the matrix factorization model. We propose a skewness-aware matrix factorization method we call SMF, which decomposes the matrix factorization framework into the basic units of rank-one matrices, and selectively combines the rank-one entries for each node pair to account for skewed network latency distributions. The first challenge is how to smartly combine the rankone matrices. Existing methods assume positive correlations between the rank-one matrices and the approximation results. We relax this assumption and consider generalized scenarios where each rank-one matrix is either positively correlated, negatively correlated or irrelevant. Then, we model these relations within a unifying combination framework, by extending the well-known orthogonal matching pursuit algorithm [41]. In each iteration, we add a new rank-one matrix into the prediction, where each entry is either added, bypassed, or decreased to account for the skewness of the network latency distribution. As most node pairs are unobserved, we need to automate the selection decision for each entry. To that end, we treat the selection decision as the rating score and the service mesh as both clients and goods, and establish the selection decision problem as a collaborative filtering task that predicts missing rating scores between a set of clients and a set of goods. Inspired by this analogy, we propose an adaptive prediction method to map the collaborating filtering scores to discrete combination choice by extending the well-known maximum margin matrix factorization method [42], [48]. Finally, extensive experiments using real-world data sets confirm that our approach reduces the relative errors by a factor of 0.2x to 10x compared to state-of-the-art methods. SMF converges fast and stably. Applying SMF to inform the low-latency detour routing [44], [45] achieves close to optimal performance. Further, we have evaluated the parameter regions where SMF obtains good approximations and verified its computational efficiency. In summary, we make three primary contributions in this paper: •We quantify the skewness of the RTT metric with respect to a set of local and global metrics, which motivates novel insights to develop skewness-adaptive matrix models. •We develop a skewness-aware matrix factorization framework SMF that automatically and adaptively selects rankone latent factors for different node pairs. •We perform extensive experiments on real-world data sets to confirm that that SMF finds a good balance between the low rankness and the adaptation to skewed RTT distributions. The rest of the paper is organized as follows. Section II summarizes the related literature on the service mesh monitoring. Next, Section III presents the background and states the requirements of adapting to skewed latent factors. Next, Section IV introduces the basic ideas of the skewness-ware matrix factorization. Section V presents the detailed algorithms and analysis of the performance and the parameter choices. Section VI reports the extensive simulation experiments compared with state-of-the-art methods. Finally, we conclude in Section VII. II. RELATED WORK Extensive studies have been made to enable network latency measurement for large-scale distributed systems and data center networks. We only introduce representative studies that are most related to us. Mesh Network Monitoring: iPlane [36] predicts end to end network latency based on an Internet topology model. iPlane issues active probes from wide-area vantage points to routable network addresses. PingMesh [25] collects allpair round-trip time (RTT) measurement system in several scales for geo-distributed data centers, which accumulates 24 TBs of probe results each day. [54] approximates network latency in an OpenFlow network environment based on control messages to and from the OpenFlow controller. Mobilyzer [39] provides a controlled and isolated library for mobile network measurement experiments. Our work is complementary to these studies by predicting missing measurements. Network Latency Prediction: In order to reduce the measurement overhead, researchers proposed to predict network latency with network coordinate methods. GNP [38] pioneers this field by modeling the network latency matrix with an Euclidean coordinate system and fitting this coordinate system via a multidimensional scaling technique. Vivaldi [11] combines a low-rank model with a height model that reflects the first-hop delay of traversing the accessing link. IDES [37] and DMFSGD [33] embed nodes using the two-factor matrix factorization. Further, DMFSGD [33] proposes a generalized framework of the low-rank matrix factorization via the SGD optimization technique. Liu et al. [34] decompose the latency metric to a distance component and a network feature component by combining the Vivaldi and the matrix completion theory. Fu et al. [19] stabilizes the matrix factorization process under churns via the relative coordinates. Zhu et al. [57] propose an adaptive matrix factorization approach by data transformation. The transformed metrics become more symmetric than the raw metrics, however, retransformed metrics become less stable when the observation is incomplete, since the estimator’s inaccuracy will be amplified exponentially with respect to the transformation base. 3 A second approach is to estimate the quantiles [59] instead of the average RTT metric based on the quantile regression framework. Generally, these methods assume that each pair of nodes shares the same set of latent factors, which may not hold when node pairs experience diverse hidden factors due to divergent routing paths, dynamic path conditions, or transient network congestions. Our work addresses this challenge via a skewness-aware matrix factorization model that selectively combines the latent factors. Matrix Completion: Our study is related with the matrix completion theory [8], which considers the problem of recovering an incomplete matrix via a subset of observed entries. For a rank-r m ×nmatrix (r(m, n)) that meets an incoherent1condition, a unique rank-rmatrix can be recovered with a high probability. When the rank is unknown, a generalized matrix completion problem seeks to best fit a partially observed matrix with a minimum rank. Unfortunately, minimizing the matrix rank exactly is NP-hard [8]. OR1MP [47] iteratively computes a rank-one matrix of the top-left and -right singular vectors of the approximation residual. However, it is generally impossible to exactly recover the SVD result for a partially observed matrix. Further, our study shows that, not all rank-one matrices are useful to reduce the approximation residual, moreover, different node pairs are correlated with heterogeneous latent factors. Traffic Matrix Interpolation: Real-world traffic matrices are usually incomplete. Consequently, interpolating missing entries becomes important. Note that traffic matrix interpolation is a related, but different problem, with different properties. Xie et al. [52], [51], [50] exploit hidden spatial and temporal structures with three-dimensional low-rank tensors, which effectively reduces the estimation error. Zhang et al. [56] interpolate incomplete traffic matrices with structure regularized low-rank matrix factorization and local interpolation procedures. LENS [10] models the traffic matrices as the sum of multiple matrices that are positively correlated with the traffic matrix. Our work is a complement to these studies by proposing a new model that keeps the low-rank interpretation and adapts well to skewed latent factors. III. PROBLEM STATEMENT In this section, we first present the measurement environment for the mesh-structured cloud services, then introduce the matrix factorization results, and discuss the open questions. A. Measurement Architecture A service mesh consists of a set of nodes located in mega data center networks or edge data-center networks. Each node hosts a set of networked micro-services, as discussed in the introduction. Service to service communication is frequent, while the latency between sending service requests and obtaining responses should meet network SLAs. As the network issues may arise in unknown locations, it is vital to maximize the network visibility between nodes for 1A well-known assumption of recovering a low-rank rmatrix is incoherence that spread the singular vectors out, which keeps the matrix not to be closely aligned with the coordinate axes. troubleshooting and proactive network optimization. To meet these needs, it is necessary to monitor network latency between any node pairs [33], [25]. Therefore, a global view of the network conditions is needed, as the service mesh is managed under a single entity. We organize the measurement system as two components: a data plane that consists of service-mesh nodes and a control plane on a logically centralized server, inspired by the software defined networks [36], [25], [54]. We assume that, the service mesh should have synchronized their clocks, as otherwise we could not correlate the network problems in different locations. The synchronization protocols such as Network Time Protocol (NTP) [2] or the IEEE 1588 Precise Time Protocol (PTP) [1] can provide millisecond-level precision for geo-distributed nodes. (i) At the control plane, the logically centralized controller schedules the RTT measurement process on the data plane. To that end, the controller randomly samples a small list of nodes as probing targets. The choices of probing targets are randomized for different nodes for load balancing. The number of probing targets depend on the measurement capability of the data plane. For a scale of hundreds of nodes, our experiments show that tens of probing targets is sufficient to obtain a good accuracy. Further, the controller handles the system dynamics, since an offline node is useless and should be detected and filtered by the controller. Accordingly, the controller keeps the online status of the data plane as volatile states in the main memory. Each online node periodically sends a heart-beating message to the centralized controller to notify its online status. After the controller successfully obtains a heart-beating message from a node, the controller piggybacks a list of sampled online nodes. The frequency of the heart-beat messages is platform-dependent. Stable platforms with dedicated nodes could choose a long period, while edge platforms should choose a relatively short period to reflect system churns. (ii) At the data plane, each service-mesh node performs a number of measurements towards other nodes in the same service mesh. It downloads the list of probing targets from the controller, and measures the round-trip time (RTT) towards these probing targets in a periodical approach. After collecting the RTT samples in an interval, each node uploads the RTT results to the persistent storage that is accessed by the controller. The data plane could use any kinds of measurement methods. For example, at the network or transport level, the data plane may choose ICMP or TCP protocol based measurement methods; while at the application level, the data plane could use RPC or HTTP protocol based methods. Generally, the RTT value amounts to the absolute difference between the time of sending a request message to the probing target and that of receiving the response message from this probing target. The unit of a measurement interval determines the granularity of the monitoring process. Increasing the sampling interval towards a probing target yields a coarser measurement granularity. Generally, the time interval should be large enough to obtain the response from the slowest probing target in the 4 system, otherwise, a data-plane node may not collect the RTT samples from some of its probing targets. B. Challenges for RTT Matrix Completion For a set of Nnodes, the pairwise RTTs between Nnodes in an interval can be represented as a N-by-Nmatrix D. The state-of-the-art approaches predict pairwise RTT values based on the matrix factorization approach, which factorizes a matrix D∈RN×Nas a product of two low-dimensional factor matrices F∈RN×rand G∈RN×r, i.e., D≈FGT, where rN, and Tdenotes the transpose of a matrix. Generally, a matrix factorization model is equivalent to a sum of a set of rank-one matrices: ˆ D=FGT=X k FkGT k(1) where Fk, Gkdenote the k-th (k≤r) column vector of the matrix Fand G, respectively. The objective function seeks to minimize the approximation residual between the observed entries and the sum of the rank-one matrices: min F,G     D− r X k=1 F∗kGT ∗k    (2) As the matrix factorization approximates the target matrix Dvia the sum of each rank-one matrix FkGT k, each of which is assumed to be positively correlated with the approximation results. Consequently, the matrix factorization assumes that all node pairs experience the same set of latent factors. Unfortunately, this assumption may not hold for the RTT metric, which additively consists of many unobservable factors, e.g., the propagation latency, the queueing latency, the transmission latency. The additive character implies that, for two node pairs with a number of different routing links, their RTTs are likely to be affected by independent latency components, or hidden factors as they are invisible from the end. Consequently, different node pairs are likely to experience a diverse set of latent factors that are unobservable from the end to end measurement, as the RTT metric only reveals the sum of all factors, not individual factors. C. Empirical Distributions Next, we empirically analyze the RTT characteristics of real-world data sets and motivate the prediction requirements. 1) Data Set Real-world network connections are heterogeneous, therefore, an ideal network latency prediction algorithm should be resilient to different network connections. We choose three kinds of publicly available data sets different in terms of scale and distributions, which provide an ideal benchmark to study the generalized performance of different prediction algorithms: •Seattle: The Seattle platform is an open peer-to-peer cloud computing platform that includes donated personal devices like personal computers, laptops and mobile phones [9] This data set was collected in summer 2014 for three hours between 99 nodes. Each interval aggregates pairwise RTT samples within 15.7 seconds, which indicates short-term dynamics between Seattle nodes. TABLE I BASIC STATISTICS OF DATA SETS. Trace Interval Mean STD Min Max Seattle 15.7s 0.37s 0.90s 0.01s 90.50s PlanetLab 14.7-hr 147.38ms 103.65ms 0.064ms 7892.8ms RIPE 2-hr 118.47ms 106.02ms 0.08ms 10,425.75ms 0 50 100 150 0 0.2 0.4 0.6 0.8 1 interquartile range (ms) Pr(X>x) (a) PlanetLab 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 interquartile range (s) Pr(X>x) (b) Seattle 0 100 200 300 0 0.2 0.4 0.6 0.8 1 interquartile range (ms) Pr(X>x) (c) RIPE Fig. 2. The CCDFs of interquartile ranges of each row vector in the dataset. •PlanetLab: The PlanetLab platform has been widely used in many previous network prediction studies [36], [33], [19]. This data set was collected in 2013 for a 9day period between 490 distributed nodes. Each interval aggregates pairwise RTT measurements within 14.7 hours, which represents long-term RTT trends between PlanetLab nodes. •RIPE: The RIPE Atlas measurement platform consists of tiny networked devices that issue measurements to a small number of addresses, most of which are chosen by the platform owner. There are two kinds of devices, i.e., Probes and Anchors, depending on their measurement capability. The Anchor is more powerful than the Probe. As we need dense RTT matrices to train the matrix factorization methods, we choose RIPE Atlas nodes that are powerful enough to probe each other. This data set2 was collected in August 19, 2018 between 250 RIPE Atlas nodes. Table I summarizes several basic statistics about the data sets. We can see that the data sets span wide ranges. 2) Distributions First, we compare the RTT distributions for each node. We calculate the interquartile ranges for each RTT vector, i.e., the difference between the 25th and the 75th percentiles of the samples in the vector. In Figure 2, we see that the interquartile ranges span wide intervals, thus different node pairs experience divergent hidden factors. Next, we evaluate the correlation between pairs of nodes. We compute the linear correlation coefficient [43] for any pair of nodes iand j, defined as c(i, j) = h~ Di◦~ Dji−h~ Diih ~ Dji σ(~ Di)σ(~ Dj), where Ddenotes the pairwise RTT matrix, ~ Djdenotes the i-th row vector, ◦denotes the hadamard operator ~ Di◦~ Djk= DikDjk for i, j, k ∈[1, N],h·i denotes the average of the vector, and σ~ Di=rD~ Di◦~ DiE−D~ DiE2 . Figure 3 plots the heat map of the linear correlation coefficients of each node pair. Darker pixels correspond to larger correlation coefficients. We can see that two to three groups 2https://data-store.ripe.net/datasets/atlas-daily-dumps/ 5 (a) PlanetLab (b) Seattle (c) RIPE Fig. 3. Heat maps of the pairwise correlation coefficients. We choose the first interval for illustration purpose. Varying the intervals yields the same conclusions. We randomly select a number of nodes from the data set as the landmarks, calculate the vector of correlation coefficients from this node to the landmarks. Then we compute the K-means clustering with these feature vectors based on the Lloyd method [53]. Finally, we reorganize the correlation coefficient matrix by putting nodes in the same cluster at adjacent positions, and plot the heat map of the matrix where darker pixels correspond to larger correlation values. During the experiments, we set the number of landmarks to 16, the number of clusters to three. The same conclusions hold as we vary the parameters. TABLE II FREQUENTLY USED NOTATIONS IN THIS PAPER. Notations Meaning NNumber of nodes D,ˆ DRTT matrix and the estimated RTT matrix Xrank-1 matrix β,ˆ βSign matrix and the estimated sign matrix EApproximation residual F,Grank-1 matrix model (~u,~v, ~ θ, b)Sign matrix model rApproximation rank SiProbing targets of node i npNumber of probing targets rsDimension of the estimated sign matrix λRegularized parameter of nodes are separable from the rest of the plot, where intragroup correlation coefficients are relatively larger than node pairs from different groups. As a result, the pairwise RTT distribution is highly skewed, and intra-group nodes are more likely to experience similar latent factors than inter-group node pairs. Accordingly, we need a powerful representative model to account for skewed latent factor distributions in the RTT matrix. IV. SKEWNESS-AWARE MATRIX FACTORIZATION Having analyzed the RTT characteristics of three data sets, the RTTs of different node pairs are likely to be affected by divergent hidden factors, which should be accounted by the prediction algorithm. Next, we first present naive approaches to account for divergent latent factors and discuss its limitations, then present a new matrix factorization model that keeps the interpretation of the low-rank representation and adapts to the skewness of the RTT distributions. Table II summarizes frequently used notations used in this paper. 0 2 4 6 8 10 0 0.1 0.2 Residual 0 2 4 6 8 10 0 0.1 0.2 0.3 0.4 Num singular values Relative error (a) PlanetLab 0 2 4 6 8 10 0 0.1 0.2 0.3 0.4 Residual 0 2 4 6 8 10 0 0.5 1 1.5 2 Num singular values Relative error (b) Seattle 0246810 0 0.1 0.2 0.3 0.4 Residual 0 2 4 6 8 10 0 0.2 0.4 0.6 Num singular values Relative error (c) RIPE Fig. 4. Energy contained in the top-ksingular values VS. the relative error of the rank-kapproximation using the SVD method. For a squared and complete matrix D, SVD calculates a Frobenius-norm optimal approximation with respect to D[24]. The SVD represents Das D=USV T, where VTdenotes the transpose of V,U∈RN×Nis an orthogonal matrix, i.e., UUT=UTU=I,V∈RN×Nis also an orthogonal matrix, i.e., V V T=VTV=I,Sis a squared diagonal matrix consisting of a vector of descending-ordered real numbers (σ1,··· , σN), where σi≥0for i≤N. A. Strawman Approaches 1) Increasing Rank A straightforward approach is to increase the rank of the matrix factorization model. In mathematics, the rank of a matrix is the maximal number of columns or rows that are linearly independent with each other, which is neatly characterized by the number of positive singular values of the SVD of this matrix [24]. The SVD can be equivalently represented as the sum of krank-1 matrices ˆ D=Pk i=1 δiuivT i, where δi is the i-th singular value, uiis the i-th left singular vector, and viis the i-th right singular vector. Each rank-1 matrix represents a latent factor that introduces a degree of freedom to approximate the matrix. We compute the relative error |Dij −ˆ Dr ij | Dij for each node pair between the rank-kapproximation ˆ Drand the RTT matrix D. Further, we compute residual fraction of the total variance captured by top-k singular values as 1−Pk iδ2 i PN iδ2 i , where δi represents the i-th largest singular value. Figure 4 shows the residual variance and the relative error with increasing number of top singular values. First, 90% of its variance can be captured with two to four top singular values, therefore, the RTT matrix is approximately low-rank. Second, the relative error is still high. A longer tail of singular values implies a higher relative error, since the residual of the SVD approximation is correlated with the remaining set of singular values, represented as D−ˆ Dk=PN i=1 δiuivT i− Pk i=1 δiuivT i=PN i=k+1 δiuivT i. 2) Incorporating Weights A second approach is to regularize the latent factor model for different node pairs. For example, let ˆ D0[i, j] = ˆ D[i, j]· w[i, j], where wdenotes the weight matrix. The estimation of each node pair is scaled independently, which introduces N2degrees of freedom to the prediction ˆ D. Accordingly, the weighted model has enough degrees of freedom to recover any matrix exactly. For example, [58] proposes to assign an Nby-Nreal-valued weighted matrix to the low-rank matrix, and estimates the weight matrix in a matrix completion framework. As [58] uses a single weight matrix to regularize every latent factor, it assumes that each latent factor is of equal importance for each node pair, which may not hold due to complex routing decisions and varying latency components. 6 Instead of assigning real-valued weights to low-rank models, we propose to decompose the matrix factorization framework into the basic units of rank-1 matrices, and selectively combine the rank-1 entries for each node pair, inspired by the wellknown orthogonal matching pursuit algorithm [41]. B. Ideal Skewness-aware Model Next, we present an ideal model to account for skewed latent factors assuming that we obtain the complete RTT matrix. In the next subsection, we relax this assumption and present a practical approach. Recall that matrix factorization is equivalent to the sum of rank-1 matrices, while each rank-1 matrix serves as a latent factor. To make the matrix factorization be aware of skewed latent factors, we should incorporate a rank-1 matrix entry for a node pair only if this node pair is correlated with this latent factor. 1) Correlation Model Let kdenote the index of the current rank-1 matrix (1≤ k≤r). Let ˆ Dkbe the approximation result of combining the first krank-1 matrices, and Ek=D−ˆ Dk−1the current residual, where ˆ D0represents an empty matrix. We define three kinds of correlation types to combine the latent factor for each node pair (i, j)as follows: (i) Positive correlation: If the approximation residual Ek[i, j]will be smaller by adding the rank-1 entry Xk[i, j]to the current approximation ˆ Dk−1[i, j], i.e., D− ˆ Dk−1[i, j] + Xk[i, j]<Ek[i, j], then we set the rank-1 entry Xk[i, j]to be positively correlated with the current approximation residual. Accordingly, we update the approximation for the node pair i, j as: ˆ Dk[i, j] = ˆ Dk−1[i, j] + Xk[i, j]. (ii) Negative correlation: If the approximation residual Ek[i, j]will be smaller if we subtract the rank-1 entry Xk[i, j]from the current approximation ˆ Dk−1[i, j], i.e., D−ˆ Dk−1[i, j]−Xk[i, j]<Ek[i, j], then we set the rank-1 entry Xk[i, j]to be negatively correlated with the approximation residual. To reduce the approximation error, we update the approximation as: ˆ Dk[i, j] = ˆ Dk−1[i, j]−Xk[i, j]. (iii) Irrelevance: Finally, if neither positive nor negative correlations hold, we should skip this rank-1 component Xk[i, j] for node pair i, j, since otherwise, the approximation residual will increase. In other words, this rank-1 entry Xk[i, j]is irrelevant with the current approximation residual. Accordingly, the approximation residual either decreases monotonically (positive and negative correlations) or keeps the current status (irrelevance). Thus we not only keep the interpretation of the matrix factorization, but also account for the skewed latent factors. 2) Algebraic Representation Model To unify the correlation model with the estimation purpose, we represent the correlation model in a compact format. As the “positive correlation, negative correlation, irrelevance” choices are equivalent to three discrete choices (+1, -1, 0), respectively, we summarize them in a pairwise matrix (denoted as a sign matrix)β∈ {+1,−1,0}N×N. Fig. 5. The diagram of the ideal correlation model. Consequently, we transform the choice of the correlation model from a combination optimization problem to a linear algebraic problem. Now we can directly represent the k-th approximation residual ˆ Ek=βk◦Xk. The overall estimation amounts to ˆ D=Pr k=1 ˆ Ek=Pr k=1 βk◦Xk. 3) Workflow The ideal model can be represented as a recurrent framework that predicts the residuals layer by layer, as illustrated in Figure 5. Each recurrence aims to find a skewness-aware rank-1 matrix to minimize the current residuals. The model sequentially finds a rank-1 matrix with respect to the residual and a sign matrix according to the correlation types. Then, it combines the rank-1 matrix and the sign matrix to approximate the residual. Finally, the overall estimation is updated and the refreshed residual is forwarded to the next recurrence. Lemma 1 shows that the ideal model either decreases the approximation residual monotonically or keeps the current status in each layer. The proof is put in the appendix. Lemma 1. Assume that we obtain the perfect sign matrix, for all k≥1,kEk+1k≤kEkk Unfortunately, manually deciding the selection choices is only possible for observed RTTs, but impossible for unobserved node pairs. To address this challenge, SMF predicts the selection choices to make the ideal model practical. C. Our Work 1) Overview SMF implements the ideal model with a layerwise learning framework. SMF unrolls the recurrent framework of the ideal model to a chained sequence of layers, and trains the model layer by layer similar to the stacking architecture of the deep neural network [27]. Each layer ktakes the partially-observed residual as the input, and computes a new rank-1 matrix by maximizing the correlation with the current residual, then combines each entry of this rank-1 matrix into the current approximation based on the sign matrix, both of which are learnt from partial observations. Next, SMF calculates the approximation residual at this layer, and finally refreshes the residual and forwards that to the next layer until reaching the final layer. SMF improves the matrix factorization in several aspects: •Explainability: SMF progressively finds a rank-1 matrix to best explain the residual in each layer. Further, higher layers’ rank-1 matrices are trained based on the residuals with respect to lower layers’ rank-1 matrices. Thus there exists no ambiguity for the rank-1 matrices. 7 •What-if Analysis: SMF optimizes a separate rank-1 matrix for each layer, and sequentially combines them with the sign matrices to produce the residual approximation. Thus SMF enables the operator to test the choice of rank based on the layerwise RTT approximation from one to r. •Robustness: SMF formulates the model with respect to the residual in each layer, which is known to be robust to the gradient-vanishing problem, as proved by the ResNet method [26]. •Modularity: SMF divides each layer to two conditionally independent optimization problems and combines layerwise estimation into a modular framework. 2) SMF Components (i) Rank-1 Matrix As each rank-1 matrix amounts to a product of two vectors [24], [47], we represent the rank-1 matrix Xkas Xk= Fk(Gk)T, where Fk, Gkdenote two vectors of length N. Moreover, to keep the latent factor interpretable, we enforce the rank-1 matrix Xkto be nonnegative by approximating the absolute-value of the residual. When the residual matrix is completely observed, the optimal vectors can be calculated based on the SVD of Ek[24], [47]: Xk=δk 1uk 1(vk 1)T, where δk 1,uk 1and vk 1denote the largest singular value, the leftmost and the rightmost singular vectors of Ek. While when the observation is incomplete, we formulate a rank-1 matrix factorization problem to estimate the vectors Fkand Gk. (ii) Sign Matrix We utilize the analogy between the sign matrix estimation and the collaborative filtering problem that recovers missing discrete rating scores (e.g., one to five stars) between a set of clients and a set of goods, and predict the sign matrix βkwith a Maximum Margin Matrix Factorization (MMMF) [42], [18]. We represent pairwise signs with a low-dimensional coordinate model. We assign each node a rs-dimensioned coordinate (~ui,~vi)and a threshold vector ~ θithat serves as boundaries to obtain the discrete signs. We make the coordinate model tolerate skewed distributions by incorporating a bias parameter to each node similar to [48], which captures the local factors in each machine. Let each node ikeep a bias variable bi, the coordinate distance Yij from node ito node jis calculated as the sum of the vector dot-product result plus the bias scalars: Yij =~ui~vj+bi+bj(3) For a pair of nodes iand j, we represent the sign from node ito node jby mapping the coordinate distance ~ui~vjto “-1”, 0, “+1” using node i’s threshold values ~ θi: •If ~ui~vj≤~ θi(1), the estimated sign b βij = -1; •If ~ θi(1) ≤~ui~vj≤~ θi(2),b βij = 0; •Otherwise, b βij = +1. 3) Example We provide a simple example to illustrate the SMF algorithm. Note that our purpose here is not to characterize the performance of the SMF. Fig. 6. Four nodes organized in a hierarchical topology. Given four nodes as shown in Figure 6 with an RTT matrix: D=0 20 80 80 20 0 80 80 80 80 0 40 80 80 40 0  For example, we build a sampled RTT matrix to model the partial choices of probing targets as follows: D0=0 20 0 80 20 0 80 0 0 80 0 40 80 0 40 0  We compute a rank-1 matrix with respect to the illustrative matrix D0, which gives F1= (−6.6791,−6.6791,−7.1096,−7.1096), G1= (−6.6791,−6.6791,−7.1096,−7.1096). The first two nodes and the last two nodes are identical with each other, respectively. Therefore, the rank-1 components captures the global proximity index. Further, F1and G1are identical, as the RTT matrix is symmetric. Multiplying F1and G1yields the rank-1 matrix: X1=0 44.6102 47.4854 47.4854 44.6102 0 47.4854 47.4854 47.4854 47.4854 0 50.5459 47.4854 47.4854 50.5459 0  Next, we calculate the residual E1with respect to the RTT matrix Das follows: E1=0−24.6102 32.5146 32.5146 −24.6102 0 32.5146 32.5146 32.5146 32.5146 0 −10.5459 32.5146 32.5146 −10.5459 0  Some RTTs are overestimated, while the others are overestimated, therefore, it is vital to selectively refine the residual. Given the partially observed residual matrix E0 1: E0 1=0−24.6102 0 32.5146 −24.6102 0 32.5146 0 0 32.5146 0 −10.5459 32.5146 0 −10.5459 0  , we find a new rank-1 matrix similar to step one, which gives: F2= (−4.7934,−4.7934,−4.3030,−4.3030),G2= (−4.7934,−4.7934,−4.3030,−4.3030). The new rank-1 matrix amounts to: X2=0 22.9771 20.6263 20.6263 22.9771 0 20.6263 20.6263 20.6263 20.6263 0 18.5161 20.6263 20.6263 18.5161 0  In order to automate the combination choices, we extract the sign matrix with respect to the partially observed residual matrix E0 1: β=0−1 0 1 −1 0 1 0 0 1 0 −1 1 0 −1 0  based on Algorithm 2, and predict the sign matrix via the MMMF framework, which gives: U:0.6074 −0.3540 0.8816 −0.1533 −1.0691 −0.2674 −0.5347 0.5162 , V :−0.8211 −0.1208 −0.6933 0.1867 0.8124 −0.9372 0.3635 0.0986 , θ :0.3265 0.3686 0.2530 0.2870 0.3905 0.5373 0.0705 0.0771  and the bias vector ( 0.1988 ,0.1186 ,-0.0130 ,0.6004). The coordinate distance matrix can be derived by Eq (3) as: ˆ Y=0−0.1698 1.0110 0.9851 −0.3880 0 0.9655 1.0243 1.0959 0.7969 0 0.1724 1.1759 1.1861 −0.3308 0  8 Fig. 7. The diagram of the SMF method. Next, we automatically infer the sign from the coordinate distance matrix via the threshold: ˆ β=0−1 1 1 −1 0 1 1 1 1 0 −1 1 1 −1 0  which accurately recovers the sign matrix. Based on the estimated rank-1 matrix and the sign matrix, we obtain the selective combination X2◦βas: X2◦β=0−22.9771 20.6263 20.6263 −22.9771 0 20.6263 20.6263 20.6263 20.6263 0 −18.5161 20.6263 20.6263 −18.5161 0  Plugging this rank-1 approximation to the rank-1 approximation, yields ˆ D2=0 21.6331 68.1117 68.1117 21.6331 0 68.1117 68.1117 68.1117 68.1117 0 32.0298 68.1117 68.1117 32.0298 0  The estimation matrix ˆ D2is much more accurate than the rank-1 estimation X1. Thus the selective combination effectively penalizes the approximation residuals. V. OPTIMIZATION METHODS Having discussed SMF’s model, we next present optimization details of the skewness-aware matrix factorization. As stated in Section III-A, the monitoring framework consists of a logically centralized controller in the control plane and distributed service-mesh nodes in the data plane. The controller predicts missing measurements in each interval: It aggregates the measurements of an interval as a partially observed matrix, and runs SMF to estimate missing matrix entries. A. Optimization Workflow We summarize the workflow of SMF in Figure 7: (i) For the first layer, let E1be the partially observed RTT matrix D, which consists of aggregated RTT samples of an interval collected from the data plane (S1). (ii) We compute a rank-1 matrix Xk(kis initialized as one for the first layer) to maximize the correlation with the residual Ek(E1is initialized to the partially observed matrix Dfor the first layer) (S2). (iii) Next, we compute a sign matrix βkto determine the combination choices of the current rank-1 matrix Xkto minimize the approximation error of the residual Ek(S3). (iv) We approximate the residual as the dot product of these two matrices (S4): ˆ Ek=Xk◦βk, and update the accumulated RTT prediction as: ˆ Dk=ˆ Dk−1+ˆ Ek=Pk l=1 Xl◦βl, where ◦represents the element-wise product operator. Next, we update the residual for each observed entry (i, j): Ek+1 [i, j] = D[i, j]−ˆ Dk[i, j]. (v) Let k←k+1. If k < r, then we move to the next layer until reaching the final layer; otherwise, we stop and output the approximation as: ˆ D=Pr l=1 Xl◦βl. Fig. 8. Graphical model of key variables in the optimization workflow. B. Loss Function Based on the optimization workflow, we see that the rank1 matrix is optimized with respect to the residual from the last layer, while the sign matrix is optimized with respect to the observed sign samples, which is calculated based on the correlation type between the residual and the approximated rank-1 matrix. Thus, each layer optimizes two conditionally independent optimization problem, as clearly shown in Figure 8. Therefore, we separate the optimization problem to decomposed subproblems and design modular optimization workflow, due to the conditional independence between the rank-1 matrix and the sign matrix. 1) Rank-1 Matrix The input to the rank-1 matrix completion is the partially observed residual. We compute the residual based on prediction ˆ Dk−1from the first to the k−1-th layer: Ek[i, j] = D[i, j]−ˆ Dk−1[i, j](4) for (i, j)∈Observed samples and k≥2. To find the rank-1 matrix for the k-th layer, we minimize the difference between the rank-1 matrix and the absolute-valued residual Ekas: min  Fk(Gk)T− |Ek| 2 F(5) We formulate a regularized objective function for Eq (5) to avoid the overfitting issue [60], [21]: min Jk= Fk(Gk)T− |Ek| 2 F+λ FT kFk−GT kGk 2 F(6) where λdenotes a regularized factor. The regularization  FT kFk−GT kGk 2 Faddresses the scaling ambiguity as Fk(Gk)T= (FkR)(GkR)Tholds for any orthogonal matrix R[60]: Suppose Jkis defined in Eq (6), any stationary point (Fk, Gk)of Jkwith ∇J(Fk, Gk) = 0 implies that FT kFk=GT kG. 2) Sign Matrix We extract the signs of observed node pairs as the input to the sign prediction problem and predict signs of for unobserved node pairs. The sign matrix is vital to keep the combined prediction monotonically decreasing. We next present a simple algorithm to set the sign matrix βthat is consistent with the selective combination rules (presented in Algorithm 2 as shown in the Appendix B): 9 First, if adding the new component reduces the absolutevalued residual, then Xij is positively correlated, so we set the corresponding sign to “+1”; second, if subtracting the new component reduces the absolute-valued residual, then Xij is negatively correlated, so we set the sign to “-1”; otherwise, we set the sign to “0”. Generally, let the ground-truth sign β[i, j]of each observed node pair (i, j) be shifted to the interval {1,2,3}for ease of performing the maximum-margin matrix factorization (MMMF). We formulate a separable soft-margin classification error LY, ~ θfor each pair of observed node pair (i,j) [42], [14], [18]: L~ui,~vj, bi, bj,~ θi= βij −1 P c=1 h~ui·~vj+bi+bj−~ θi(c)+ 2 P c=βij h~ θi(c)−(~ui·~vj+bi+bj) = 2 P c=1 hTc ij [βij]·(θic −(~ui·~vj+bi+bj))(7) where θic =~ θi(c),Tc ij [βij] = n+1 c≥βij −1c<βij serves as an indicator function for the sign value βij , and h(z) = 1 2−z z<0 0z>1 1 2(1−z)2otheriwse represent a smoothed hinge loss with derivative: h0(z) = n−1z<0 0z>1 z−1otheriwse Thus, we can optimize the loss function with gradient based optimization techniques. Further, We present a regularized loss function with respect to Eq (7) for each node i[42], [14], [18] for robustness against outliers and missing items: gk(xi) = X ∈Si L~ui,~vj, bi, bj,~ θi+λk~uik2 F+k~vik2 F+kbik2 F (8) where xi=h~ui;~vi;bi;~ θii,λdenotes the regularization constant. C. Algorithm The key building block of each layer is the optimization of the rank-1 matrix and that of the sign matrix, both of which belong to the non-convex matrix factorization problem. 1) Rank-1 Matrix Optimization First, we seek a new rank-1 matrix Xk, which needs to be maximally correlated with the approximation residual. We compute the rank-1 matrix with the stochastic gradient descendent method that onverges fast to the vicinity of the optimum regions and helps escape the local minimum for distributed optimization [22], [31], [46], [29], [23], [21]. We optimize Eq (6) with the stochastic gradient descendent method. Let Fk(i)and Gk(i)denote the row vectors of node i. We optimize each node i’s row vectors in Trounds, whereas in each round t, we compute the partial gradients of the row vectors with respect to probed targets of node i, and then adjust the row vectors towards the approximated negative gradient direction scaled by a learning-rate parameter η1, which is automatically chosen using the line search method [33]. We summarize the optimization steps in Algorithm 3 (Appendix B). Algorithm 1: Optimization procedures of the SMF algorithm. 1SMF(D,r) input : Partially observed RTT matrix D, rank r output: Estimated RTT matrix ˆ D 2Compute a rank-1 matrix X1=F1(G1)Tthat is maximally correlated with the partially observed RTT matrix Dby calling RankOne(D)in Algorithm 3 (Appendix B); 3Set the RTT prediction as ˆ D1=X1; 4for Integer k∈[2, r]do 5Compute the partially observed residual Ek=D−ˆ Dk−1; 6Find a rank-1 matrix Xk=Fk(Gk)Twith respect to Ek by calling RankOne(|Ek|)in Algorithm 3 (Appendix B); 7Calculate the partially observed sign matrix βkby calling SignRule(Ek[i, j],Xij )in Algorithm 2 (Appendix B); 8Find a sign matrix ˆ βkby calling SignMatrix(βk)in Algorithm 4 (Appendix B); 9Update the RTT prediction ˆ Dk=ˆ Dk−1+Xk◦ˆ βk; 10 Set k←k+ 1; 11 return ˆ D= r P k=1 Fk(Gk)T◦ˆ βk; 2) Sign Matrix Optimization Second, we seek a sign matrix that determines the correlation types for each node pair between the residual and the current rank-1 matrix. We compute the sign matrix based on the MMMF framework [42] that adaptively penalizes the difference between the estimated signs and the ground-truth ones. We optimize Eq (8) with the nonlinear conjugate-gradient optimization [42], [18]. We form a vector xias the concatenation of the coordinate components of i, i.e., xi(s) = h~ui;~vi;~ θi;bii, and adjust each node’s vector xiin Trounds, whereas in each round s(s≥1), we compute the conjugate direction of the vector xiwith respect to probed targets, and move the vector xia small distance in the direction of the conjugate direction scaled by a learning-rate parameter η2, which is automatically computed to meet the Strong wolfe line search conditions [12]. The optimization steps are summarized in Algorithm 4 (Appendix B). 3) Putting It All Algorithm 1 summarizes the overall workflow. First, we find a rank-1 matrix that is maximally correlated with the partiallyobserved RTT matrix D(In line two). Next, we set the current approximation based on the rank-1 matrix X1( in line three). Then, we perform layerwise training (from lines four to ten), by finding the rank-1 matrix and the sign matrix, and forward the residual to the next loop cycle. D. Analysis 1) Convergence SMF generalizes the truncated SVD that recursively finds an approximately optimal rank-1 matrix to the residual. The SVD sums up the rank-1 matrix to form the low-rank prediction, while SMF flexibly combines the rank-1 matrix to be aware of the skewness of node pairs. As a result, the truncated SVD is a special case for SMF. 16 pair, the expected number of correctly predicted signs amounts to r(1 −), the variance is r (1 −), and for t > 0, the number of correct estimations is not within tpr(1 −)from its expectation is at most  t2. Proof. Assume that the sign matrix can be correctly predicted for at least (1 −)(for ∈[0,1]) fractions of entries in each iteration, and that the incorrect signs are spread uniformly at random over the matrix. Then the sign of any server pair is incorrectly predicted with a probability . For any matrix entry (i, j)and the index k,kEk+1 [i, j]k ≤ kEk[i, j]kholds with a probability (1 −). Since the sign matrices of different iterations are independent with each other, the always-failure probability is at most r, which decreases exponentially to zero with increasing rank r. For example, if = 0.1and r= 8, the always-failure probability amounts to r= 0.18≈0. Generally, let Zk[i, j]be an indicator function of whether the estimated sign is correct for the k-th iteration: Zk[i, j] = 1ˆ β[i, j] == β[i, j] 0ˆ β[i, j]6=β[i, j] We can see that the probability Pr(Zk== 1) = 1 −. Let Z[i, j]be the random variable of the number of correct predictions for the server pair (i, j). Then we have Z[i, j] = r X k=1 Zk[i, j](19) The expected number of correct predictions for the server pair (i, j)can be calculated as: E[Z[i, j]] = E[Pr k=1 Zk[i, j]] =Pr k=1 E[Zk[i, j]] =Pr k=1 (1 ×(1 −)) =r(1 −) (20) Further, the variance of the number of correct predictions can be written as: V ar (Z[i, j]) = V ar [Pr k=1 Zk[i, j]] =Pr k=1 V ar [Zk[i, j]] =Pr k=1 Eh(Zk[i, j])2i−(E[Zk[i, j]])2 =Pr k=1 12×(1 −)−(1 −)2 =r (1 −) (21) By the standard Chebyshev’s inequality, for δ > 0, we have P(|Z[i, j]−r(1 −)|> δ)≤r (1 −) δ2 By taking δ=tpr(1 −), where t > 0, we have P|Z[i, j]−r(1 −)|< tpr(1 −)≥1− t2(22) which implies that the probability that a particular Z[i, j]not within tpr(1 −)from its expectation is at most  t2.