A De-singularity Subgradient Approach for the Extended Weber Location Problem

1 Introduction

The extended Weber location problem is a classical optimization problem [Ostresh, 1978; Brimberg and Love, 1993; Beck and Sabach, 2015] that has been introduced to some new machine learning scenarios recently [Aftab et al., 2015; Li et al., 2016; Lai et al., 2018c; Lai et al., 2020; Lai and Yang, 2023]. It finds the point that minimizes the q-th power of the Euclidean distances from m fixed data points x1, · · · , xm ∈ Rd [Weber, 1909; Cooper, 1968;

* Corresponding author.

** The supplementary material and codes for this paper are available at https://github.com/laizhr/qPWAWS.

Chen, 1984]:

where ∥ · ∥q denotes the q-th power of the Euclidean distance, and Cq(y) is the q-th power cost function at the point y ∈ Rd . The q-th power median x∗ (especially when 1 ⩽ q < 2) has recently been found useful in rotation averaging and gives more robust results when outliers exist [Aftab et al., 2015]. We will further show its advantage in online portfolio selection [Li et al., 2016; Lai and Yang, 2023] in this paper.

1.1 The Singularity Problem

In order to have a quick insight into the aim of this paper, we first compute the gradient of Cq(y) w.r.t. y:

It is well-defined when q ⩾ 2 and one could easily find some gradient descent methods to work out the minimum [Afsari et al., 2013], which is beyond the main concern of this paper. However, it is singular at the data points {xi}mi=1 when 1 ⩽ q < 2, since at least one of the weights ∥y-xi∥q-2 becomes undefined (see Figure 1). This paper mainly addresses this case and we assume 1 ⩽ q < 2 in the rest of this paper if not specified.

Since Cq(y) is convex and continuous in Rd including these data points {xi}mi=1 , we can still turn to the subgradient approach for solutions.

Definition 1 ([Rockafellar and Wets, 2009]). The Fr´echet subdifferential of Cq at y, denoted by ∂Cq(y), is the set of all vectors v ∈ Rd satisfying

Since Cq is also proper (dom (Cq) ̸= ∅ and Cq(y) ̸= -∞, ∀y ∈ Rd ), its Fr´echet subdifferential is equivalent to the classical subdifferential in the literature of convex analysis.

Definition 2 (Subdifferential in Convex Analysis).

If Cq is differentiable at y, ∂Cq(y) reduce to a gradient ∇Cq(y). Based on the properties of Cq (1 ⩽ q < 2), it is easy to verify that at least one Fr´echet subgradient exists at each singular point xk ∈ {xi}mi=1 by satisfying the limit inferior inequality in (3). However, few existing works have been done to deal with this singularity problem, because it is usually considered a rare event in practice or just overlooked.

One of the most typical and widely-used gradient approaches is the general q-th Power Weiszfeld Algorithm (qPWA, [Cooper, 1968; Chen, 1984; Aftab et al., 2015]). Since Cq(y) is convex, y is the minimum of Cq(y) if and only if ∇Cq(y) = 0. It leads to the following update formula of qPWA at the p-th iteration:

where y(p) is the current iterate and y(p+1) is the next iterate computed by a re-weighted sum of the data points {xi}mi=1 . The weights of {xi}mi=1 are determined by the current distances fromy(p) to {xi}mi=1 and normalized by their sum. (5) can be further transformed to:

where ϑ = 1/(q ∑mi=1 ∥y(p) - xi∥q-2 ) is an automatic step size for the gradient descent approach. Therefore, the Weiszfeld algorithm (5) is actually a gradient descent method (6), which shall fail if y(p) hits one of {xi}mi=1 and 1 ⩽ q < 2.

1.2 The Significance and Difficulty of This Singularity Problem

This singularity problem actually happens frequently and unexpectedly. We can easily give a simple example. Let q = 1.1

and 6 data points be

The starting point y(0) ≜ (1.68645, 0) is chosen distinct from {xi}6i=1 . The existing algorithm qPWA gets stuck with just 1 iteration since it hits one data point x3 = (1, 0) and cannot proceed. The resulted cost function is Cq = 9.4201, which is not the true minimum. On the contrary, the proposed algorithm qPWAWS successfully escapes from the singular point x3 and finds the true solution (0, 0) with the true minimum Cq = 8.2871. [Chandrasekaran and Tamir, 1989; Vardi and Zhang, 2000] further point out that the bad starting point set {y(0)} itself may constitute a continuum set that can be dense in an open region of Rd with d ⩾ 2, even the entire Rd . Hence this singularity problem is quite significant, if not serious.

This singularity problem cannot be radically eliminated by straightforward treatments, such as perturbations and random restarts. A treatment is to find a starting point with a smaller cost than the costs of any data points {xi}mi=1 by computing {Cq(xi)}mi=1 at first [Aftab et al., 2015], but it is rather exhaustive especially when the number of data points m is large. It requires at least m iterates to find a starting point, but sometimes we only need fewer iterates to find a minimum point. Table 3 shows that our method requires only 9.31 iterates in average to find the minimum in the "NYSE(N), m = 10, q = 1.9" case, while the above treatment conducts m = 10 iterates only to find a starting point.Furthermore, if xk is exactly the minimum, we have to check whether 0 ∈ ∂Cq(xk). After all, once we figure out ∂Cq(xk), we can get rid of these troubles. In fact, there is no need to introduce additional treatments because our de-singularity approach can immediately replace the ordinary gradient step without increasing computational complexity, which is a great advantage.

In the q = 1 case, there is a remedy that removes the singular term in the gradient [Kuhn, 1973; Vardi and Zhang, 2000], but it does not investigate the subgradient ∂Cq(xk). More importantly, there are some incomplete statements in the proof of convergence, which will be fixed in Section 3.4. In the 1 < q < 2 case, neither has the subgradient ∂Cq(xk) been figured out nor has the proof of convergence been completed, yet.

1.3 Our Results

We mainly establish a complete "de-singularity" subgradient approach for the extended Weber location problem (1) (1 ⩽ q < 2). The key contributions fall into four aspects:

1. By removing the singular term, we propose a desingularity subgradient of Cq at each singular point xk ∈ {Cq(xi)}mi=1 . It can replace the ordinary gradient without increasing computational complexity.
2. If xk is not a minimum point, we pull the stuck iterate in the de-singularity subgradient descent direction that can reduce the cost Cq, so that the iterates can monotonically converge to the exact minimum. In the rest of this pa-

per, we call the new algorithm q-th Power Weiszfeld Algorithm without Singularity (qPWAWS).

1. We present a complete proof of convergence for qPWAWS, which has fixed some incomplete statements of the proofs for some previous Weiszfeld algorithms.
2. If one of {xi}mi=1 is theminimum, we prove that qPWAWS enjoys a superlinear convergence for the iteration sequence when 1 < q < 2, which is a new theoretical result. If none of {xi}mi=1 is the minimum, we show that qPWAWS enjoys a reasonable rate of linear convergence for the iteration sequence with computational experiments.

The proposed algorithm qPWAWS can use any starting point in Rd and guarantee a monotonic convergence to the exact minimum. Moreover, this de-singularity subgradient approach can also be adopted in other gradient-related algorithms for problem (1).

2 RelatedWorks on theWeiszfeld Algorithms

In this section, we assume that the data points {xi}mi=1 are distinct and non-collinear (i.e., they do not lie in a hyperplane). Thus the cost function Cq(y) in (1) is strictly convex and has a unique minimum. We review some related works on how to derive the Weiszfeld algorithms.

2.1 Convergence in the Non-singularity Case

For the basic qPWA (5), it has been proven that Cq(y(p+1)) < Cq(y(p)) if y(p+1) ̸= y(p), and that as long as the sequence of iterates does not hit any data points, it monotonically converges to the exact minimum [Kuhn, 1973; Cooper, 1968; Chen, 1984; Aftab et al., 2015], while some incomplete statements are to be fixed in Section 3.4. The rate of convergence for qPWA (1 < q < 2) has not been deduced, yet.

2.2 Iterative Re-weighted Least Squares (IRLS) Interpretation

Based on the current iterate y(p), a weighted 2-nd power cost function ˜ Cq(y) can be created to approximate Cq(y):

By setting the gradient of ˜ Cq(y) as zero, we have the same update formula as (5). It is a kind of Iterative Re-weighted Least Squares (IRLS) techniques [Chartrand and Yin, 2008; Daubechies et al., 2010; Eldar and Mishali, 2009; Aftab et al., 2015].

2.3 1-st PowerWeiszfeld Algorithmwithout Singularity

A remedy for the L1-median is studied in [Kuhn, 1973; Vardi and Zhang, 2000; Ostresh, 1978]. First, we define a de-singularity Weiszfeld transform that excludes the singular point as follows:

Then the next iterate is the combination of ˜T(y(p)) and the current iterate y(p):

where 0 ⩽ λ ⩽ 1 is a mixing parameter.

In order to drag the iterate out of the singular point and reduce the cost simultaneously, we adopt the 1-st power desingularity subgradient (a formal definition will be given in Definition 5) and set λ as follows:

Compared with (2), ∇D1(y(p)) removes the singular term in ∇C1(y(p)). If y(p) does not hit any data point, then (9) is just the same as (5). Else, λ ensures C1(y(p+1)) ⩽ C1(y(p)) and the iterates monotonically converge to the exact minimum, which are explained by [Vardi and Zhang, 2000].

3 q-th PowerWeiszfeld Algorithmwithout Singularity

Although there is a remedy for the 1-st power case, it is nontrivial to deal with the general q-th power (1 < q < 2) case. We will see that the 1-st power case and the q-th power (1 < q < 2) case are different in the characterization of minimum, the minimizing strategy, and the rate of convergence, since they have different smoothness. Nevertheless, we still present a unified q-th power (1 ⩽ q < 2) framework for the integrity of the approach.

In the rest of this paper, we generalize the q-th power (1 ⩽ q < 2) cost function in (1) to:

where ξi > 0 is the multiplicity of the data point xi. It can deal with duplication and collinearity of data points. For example, if x1 = x2 = x3 and the other points are distinct, then we can merge x1, x2, x3 and set ξ1 = 3, ξi = 1(i ⩾ 4). The collinear data points can also be merged to be non-collinear ones by the same way. Hence, we assume that the data points are distinct and non-collinear in the rest of this paper. Thus Cq(y) is strictly convex on y and there is a unique minimum

point M. To be convenient for illustrations, we set ηi ≜ ξ 1 q i and further change (12) to:

The illustration of qPWAWS consists of 5 steps:

1. The general update formula is introduced as an easy beginning.
2. The q-th power de-singularity subgradient is defined to characterize the subgradients and the minimum.

1. The q-th power de-singularity subgradient is adopted to get out of the singular point and reduce the cost simultaneously.
2. The proof of convergence is conducted.
3. The theoretical rate of convergence for the iteration sequence in a special case is deduced with the properties of the q-th power de-singularity subgradient.

3.1 General Update Formula without Coincidence

First, we consider the simplest case where the current iterate y(p) does not coincide with the data points {xi}mi=1 . We first state the general iterative update formula in this case:

As a more general case than [Kuhn, 1973; Cooper, 1968; Chen, 1984; Aftab et al., 2015], we present the following non-increasing theorem:

The proof is provided in Supplementary B.1.

3.2 Characterization of Subgradients and Minimum

Second, we establish the de-singularity subgradient and characterize the minimum point M of Cq(y) in (13). The following corollary is a direct result of Theorem 3.

Corollary 4. If y(p) /∈ {xi}mi=1 , then T1(y(p)) = y(p) ⇔ y(p) is the minimum point M of Cq(y).

The proof is provided in Supplementary B.2. If y(p) ∈ {xi}mi=1 , then the minimum characterization relies on the subgradient(s) in ∂Cq(y(p)). Without loss of generality, suppose y(p) = xk.

Definition 5 (q-th Power De-singularity Subgradient). Let Dq(y) be the main component of Cq(y) that excludes the term ηqk ∥y - xk∥q . Then the q-th power de-singularity subgradient ∇Dq(y) of Cq(y) is:

Compared with the ordinary gradient (2), the de-singularity subgradient (16) does not require more computation but checking whether y = xk for the k-th summand, which can be done at the same time and in the same loop with summand computing. Hence the de-singularity subgradient will not increase computational complexity.

Theorem 6 (Characterization of Subgradients and Minimum). Let y(p) = xk. Then

According to Fermat's rule, xk is the minimum point M if and only if 0 ∈ ∂Cq(xk).

The proof is provided in Supplementary B.3. It is obvious that ∇Dq(xk) ∈ ∂Cq(xk) for all 1 ⩽ q < 2. This is why ∇Dq(xk) is termed a de-singularity subgradient: it removes the singularity and then becomes a subgradient. One can also see that (10) is a special case of (17) with y(p) = xk and u = 0.

-∇Dq(xk) can be interpreted as the resultant implemented on xk towards other data points (see Figure 2). When q = 1 and ∥∇D1(xk)∥ is smaller than the intrinsic force ηk implemented on xk, the system remains still and Cq(xk) has reached its minimum. But when 1 < q < 2, ∥∇Dq(xk)∥ should be zero to keep the system still, no matter how large the intrinsic force ηk is.

We have characterized the minimum no matter whether the current iterate y(p) hits the data points {xi}mi=1 or not. If y(p) has not reached the minimum point M, then another iteration should be implemented. When y(p) /∈ {xi}mi=1 , the general formula (14) can be adopted. When y(p) ∈ {xi}mi=1 , we will show later that y(p) - λ∇Dq(y(p)) can get away from the singular point and reduce the q-th power cost simultaneously.

3.3 Update Formula with Coincidence

The update formula with coincidence y(p) = xk should be set up separately for q = 1 and 1 < q < 2. The q = 1 case is provided in Supplementary A.1. Now we turn to solve the 1 < q < 2 case. The following theorem shows a way to reduce the q-th power cost when the iterate gets stuck in the singular point.

Theorem7 (De-singularity Subgradient Descent Method). If y(p) = xk and ∥∇Dq(xk)∥ > 0, then there exists a λ∗ > 0 such that for any 0 < λ ⩽ λ∗, Cq(xk - λ∇Dq(xk)) < Cq(xk), 1 < q < 2.

The proof is provided in Supplementary B.4. The key point of Theorem 7 is that the negative de-singularity subgradient part (resultant) -λ∥∇Dq(xk)∥2 is of lower order infinitesimal than the positive singular part (resistance) ηqk λq ∥∇Dq(xk)∥q . Thus when λ is sufficiently small, the resultant overcomes the resistance and drags the current iterate y(p) out of the singular point and the q-th power cost is reduced simultaneously. An algorithm can be designed based on this mechanism. By omitting λ2 G(xk) + o(λ) on the left side of (43), we can approximately start with the following λ0:

As long as Cq(xk - λw∇Dq(xk)) ⩾ Cq(xk), we reduce λw with a factor ρ < 1: λw+1 ← ρλw, w = 0, 1, · · · , until we find some λ∗ such that Cq(xk - λ∗∇Dq(xk)) < Cq(xk) (see Figure 3). According to Theorem 7, this λ∗ is sure to be found. Then the next iterate is:

This strategy is very efficient to get away from the singular point: it takes only 2.75 iterates in average in the worst case of our experiments (see Table 1). It can be even more efficient if we further reduce λ0 and ρ. Besides, the choice of λ∗ can also be absorbed in momentum acceleration methods, such as Nesterov [Nesterov, 1983; Sutskever et al., 2013] and Adam [Kingma and Ba, 2015]. Details can be found in the code link provided in the section of acknowledgments.

3.4 AComplete Proof of Convergence

From (14) and (19), the q-th Power Weiszfeld Algorithm without Singularity (qPWAWS) can be designed as (1 < q < 2):

where T1(y(p)) is the update formula without coincidence (14) and T2(y(p)) is the update formula with coincidence (19). Note that the q = 1 algorithm without singularity has been proposed by [Kuhn, 1973; Vardi and Zhang, 2000], but their convergence proofs are incomplete. Thus we focus on the proof of the 1 < q < 2 case as well as fix their incomplete statements in this subsection. In calculation, the q = 1 case and the 1 < q < 2 case can be unified.

Starting with any initial point y(0), the qPWAWS generates a sequence of iterates y(0), y(1), · · · , y(p) · · · . It is necessary to prove that the sequence converges to the minimum point M unless some y(p) hits M.

Lemma 8. The sequence {y(p)} generated by qPWAWS (20) visits each xk ̸=Mat most once and will not get stuck. Except for at most a finite set of iterates, {y(p)} and M lie in the convex hull of the data points {xi}mi=1 .

The proof is provided in Supplementary B.5. It is also possible that the sequence {y(p)} converges to the data points {xi}mi=1 . In fact, T1(y) is continuous in almost the whole Rd except for the data points {xi}mi=1 where it is discontinuous. However, {xi}mi=1 are removable discontinuities. Less general versions of this characteristic are mentioned in [Kuhn, 1973; Aftab et al., 2015] but without proof. Thus we present the following lemma:

Lemma 9.

The proof is provided in Supplementary B.6. Based on this lemma, the operator T1 can be extended to the removable discontinuous data points:

When y gets into some neighborhood of a data point xk ̸= M, the operator T1 will eventually drive it out of the neighborhood. This property of T1 helps to find the real minimum point. The following lemma is a nontrivial extension of a similar conclusion for the 1-st power case in [Kuhn, 1973].

Lemma 10. If xk ̸= M, then there exists some δ0 > 0 such that for all y ∈ B(xk, δ0), Ts-1 1 (y) ∈ B(xk, δ0) and Ts1 (y) /∈ B(xk, δ0) for some s. B(xk, δ0) denotes an open δ0-ball centered at xk. Ts1 means implementing T1 for s times, and s depends on y.

The proof is provided in Supplementary B.7.

Theorem 11 (Convergence Theorem). Starting from any initial point y(0), if the sequence {y(p)} generated by qPWAWS (20) does not hit M, then y(p) → M. If y(p) hits M, the characterization of minimum (Theorem 6) ensures that this hit could be recognized and the algorithm would be stopped.

The proof is provided in Supplementary B.8.

Some notes: 1. In [Kuhn, 1973; Vardi and Zhang, 2000], the proof of convergence assumes that

limv→∞y(pv) = xk but T1(y(pv)) /∈ B(xk, δ0) for all v. However, Lemma 10 indicates that T1 may be implemented for s > 1 times to drive y(pv) out of B(xk, δ0). Besides, this s also depends on y(pv) and the uniform upper bound of s for all the y(pv) may not exist. Thus it may not be concluded that limv→∞ ∥T1(y(pv))-xk∥ ∥y(pv)-xk∥ = ∞ and it may fail the conclusions in [Kuhn, 1973; Vardi and Zhang, 2000]. Moreover, [Kuhn, 1973; Vardi and Zhang, 2000] have only shown that a subsequence {y(pv)} converges to M, but not the whole sequence {y(p)}. The convergence of the whole sequence {y(p)} should depend on the uniqueness of M.

2. In [Aftab et al., 2015], the proof of convergence does not include Lemma 10. Thus it may not guarantee that a subsequence y(pu) ∈ B(xk, δ0) and T1(y(pu)) /∈ B(xk, δ0). In this case, (57), (58), (59) and (60) may not be deduced and the contradiction may not exist.

Therefore, the proof of convergence in this paper has also fixed some incomplete statements of some related works before.

3.5 Rate of Convergence

It is also important to analyze the rate of convergence for qPWAWS when y(p) → M. With the proposed desingularity subgradient ∇Dq in (16), we can deduce the exact rate of convergence for qPWAWS in the special case M ∈ {xi}mi=1 . To the best of our knowledge, this is a new theoretical result that qPWAWS enjoys a superlinear convergence for the iteration sequence in this case. Without loss of generality, we assume M = xk. Since the subgradients are different between the q = 1 case and the 1 < q < 2 case (Theorem 6), we further divide the analysis into two parts. The q = 1 case has been deduced by [Ostresh, 1978], while the 1 < q < 2 case is somewhat complicated. The key technique for the 1 < q < 2 case is to check the order of infinitesimal of ∥∇Dq(y(p))∥ with y(p) → xk at first.

Lemma 12. If 1 < q < 2, the order of infinitesimal of ∥∇Dq(y(p))∥ is no higher than that of ∥y(p) - xk∥ with y(p) → xk. Precisely, there exists some ζ ⩾ 0 such that

The proof is provided in Supplementary B.9.

Theorem13. IfM= xk for some k, the rate of convergence for qPWAWS is:

It is a superlinear convergence when 1 < q < 2 and no worse than a sublinear convergence when q = 1.

The proof is provided in Supplementary B.10. As for the general case M /∈ {xi}mi=1 , the rate of convergence is unknown at present according to our knowledge. It is difficult to eliminate the infinitesimal ∥y(p) -M∥ and obtain a uniform upper bound of limy(p)→M ∥y(p+1)-M∥ ∥y(p)-M∥ that is no greater than 1. When q = 1, [Ostresh, 1978] shows some computational evidence that the rate of convergence for M /∈ {xi}mi=1 is usually somewhat less than that of M ∈ {xi}mi=1 . We will give some computational results for the 1 < q < 2,M /∈ {xi}mi=1 case in Section 4.3 to show that qPWAWS enjoys a reasonable rate of linear convergence. We summarize the whole qPWAWS in Supplementary A.2, which gives a complete procedure to deal with different situations that may occur in real applications.

4 Experimental Results

Since 1PWA and qPWA have already been verified to be effective in many applications [Hartley et al., 2011; Hartley et al., 2013; Fletcher et al., 2009; Yang, 2010; Huang et al., 2016], we do not intend to repeat their experiments. Instead, we focus on validating that qPWAWS successfully solves the singular problem and enjoys a reasonable rate of linear convergence.

We conduct experiments on an interesting machine learning application: online portfolio selection (OPS, [Li et al.,

2016; Lai et al., 2018c; Lai et al., 2020; Lai and Yang, 2023]). By using the notations in this paper, a data point xi ∈ Rd indicates a price vector of d assets on the i-th day. {xi}mi=1 contains the asset prices on the most recent m days. Following [Huang et al., 2016], we fed the q-th power median

to a portfolio optimization model [Huang et al., 2016] to decide the future portfolio. We adopt the NYSE(N) data set [Li et al., 2013] and propose a new CSI300 data set.

• NYSE(N): it contains daily price relative sequences of d = 23 stocks from New York Stock Exchange during 1/Jan/1985∼30/Jun/2010 (T = 6431 days).
• CSI300: it contains daily price relative sequences of d = 47 stocks from the CSI300 constituents 1 of Shanghai Stock Exchange and Shenzhen Stock Exchange in China during 16/Mar/2015∼19/May/2017 (T = 534 days).

Experiments consist of four parts: first, we verify that qPWAWS solves the singularity problem and obtains the same results as qPWA. Second, we measure the computational cost and the number of iterations for convergence of qPWAWS. Third, we give computational results on the rate of convergence for qPWAWS. Fourth, we show that the q-th power median (1 < q < 2) is more effective and robust than the 1-st power median and the 2-nd power median in some cases, thus the q-th power median is useful and it is important to solve the singularity problem. As for the parameters, if not specified, the observation window size is m = 5, which is consistent with previous related methods [Huang et al., 2016; Lai et al., 2018a; Lai et al., 2018b; Lai et al., 2022]. The tolerance threshold is Tol = 10-9 . The reducing factor is ρ = 0.1, which is a moderate value. As the observation window moves from t = 1 to t = T-m+1, there are a total number of (T -m+1) sets of data points {xi}mi=1 . Thus the q-th power median computation can be conducted for (T -m+1) times to evaluate the average performance of qPWAWS.

4.1 Addressing the Singularity Problem

If the sequence of iterates does not hit the data points, then qPWAWS is equivalent to qPWA. Without loss of generality, we examine a specific case: let the starting point y(0) = x1 for qPWAWS and y(0) = 1m ∑mi=1 xi for qPWA. Besides, we also watch to avoid getting stuck for qPWA, since it cannot deal with singularity.

To evaluate the efficiency of qPWAWS, we measure how many iterates that qPWAWS takes to successfully get away from the singular point with a reduced cost (implementing only Step 11 in Algorithm 1). For each data set, we use two different sizes of observation windows m = 5 and m = 10, and let q = 1.1 ∼ 1.9. The mean and the standard deviation (STD) of the number of iterates are shown in Table 1. qPWAWS successfully gets away from the singular points in all the experiments. As q increases from 1.1 to 1.9, the average number of iterates decreases from about 2.54 to 2 for m = 5 and from about 2.75 to 2 for m = 10 on CSI300. In fact, a smaller ρ and a smaller λ0 in (18) can further reduce the number of iterates and Step 11 only need to compute the q-th power cost. In general, the number of iterates to get out of the singular point is small.

1 http://www.csindex.com.cn

Next, we need to verify that qPWAWS obtains the same q-th power median as qPWA if the latter does not get stuck. The experimental procedure is similar to the above. In each q-th power median computation, denote the q-th power medians of qPWAWS and qPWA by ˆxWAWS and ˆxWA, respectively. Then we compute the maximum relative difference ∥ˆxWAWS - ˆxWA∥/∥ˆxWA∥ of all the (T - m + 1) times of computations, shown in Table 2. Since the differences are close to zero (< 1e - 07), qPWAWS obtains nearly the same q-th power median as qPWA.

4.2 Computational Costs and Convergence

A computer with an Intel Core i7-6700 CPU and a 4GB DDR3 memory card is used to record the computational time of qPWAWS. (T -m+1) times of computations are recorded and the mean time cost for each q and m is shown in Table 3. All the computations successfully converge and the average number of iterates for convergence (including the escaping iterates) is also shown in Table 3. All the time costs (< 1e - 03s) and the numbers of iterates (< 27) are small, which indicates that qPWAWS runs fast.

4.3 Rate of Convergence (Computational)

In Section 3.5, we have deduced the theoretical rate of convergence for qPWAWS in the special case M ∈ {xi}mi=1 . Now we conduct computational experiments to show that qPWAWS enjoys a reasonable rate of linear convergence in general. The procedure is the same as the above subsections and (T -m+1) times of computations are implemented. In each computation, we take ∥y(o-1)-y(o)∥ ∥y(o-2)-y(o)∥ as an approximation for the rate of convergence, where y(o) denotes the last iterate (i.e., the minimum) in Algorithm 1. The means and the STDs of the rates of convergence for qPWAWS with different values of q are shown in Table 4. When q changes from 1 to 1.9, the empirical rate of convergence decreases from about 0.36 to 0.06, thus the newly-developed 1 < q < 2 algorithm enjoys better rate of convergence than the long-used q = 1 algorithm. Besides, all the situations with 1 ⩽ q < 2 achieve reasonable rates of linear convergence.

In order to see how the rate of convergence changes as qPWAWS runs, we plot the sequences of rates { ∥y(p+1)-y(o)∥ ∥y(p)-y(o)∥ }o-2 p=1 with q = 1.1 ∼ 1.9 (one instance for each q) and m = 5 in Figure 4. Each plot has a slowlydecreasing period and drops sharply in the last few iterates, thus the inflection point is close to the rate of convergence for qPWAWS, which suggests a linear convergence. Besides, the plots with smaller qs are above the plots with larger qs, which accords with the results in Table 4. Other instances of rates of convergence show similar patterns and we need not plot all of them.

4.4 Investing Performance

We apply the q-th power median and qPWAWS to the OPS problem as a price prediction strategy and evaluate the investing performance on the two data sets NYSE(N) and CSI300. Interested readers are referred to [Huang et al., 2016] for the specific model setting. We change q from 1 to 2 to see how it affects the investing performance. We set m = 5 to be consistent with the window size of previous works.

The first indicator of evaluation is the final cumulative wealth (CW) when a strategy goes through the whole investment. The second indicator is the daily Sharpe Ratio (SR,

[Sharpe, 1966]), which is a risk-adjusted average return that considers both return and risk in the investment. The results in Table 5 indicate that q = 1.3 and q = 1.6 achieve the best CW and SR on NYSE(N) and CSI300, respectively. They improve on the trivial q = 1 and q = 2 to some extent, which suggests that the q-th power median (1 < q < 2) is more effective and robust than the 1-st power median and the 2-nd power median in some cases. Hence the q-th power median and qPWAWS are useful in this machine learning scenario.

5 Conclusions and FutureWorks

This paper mainly establishes a novel de-singularity subgradient approach and a corresponding algorithm (q-th Power Weiszfeld Algorithm without Singularity, qPWAWS) for the extended Weber location problem. We characterize the subgradient(s) and the optimality of any given singular point. If this singular point is not optimal, the algorithm can escape from it and reduce the cost simultaneously. This advantage makes the sequence of qPWAWS monotonically converge to the exact minimum. A complete proof of convergence for qPWAWS is also presented, which has fixed some incomplete statements of the proofs for some previous Weiszfeld algorithms. qPWAWS enjoys a superlinear convergence for the iteration sequence in the special case where the minimum point is a singular point, which is a new theoretical result.

Experiments with real-world financial data sets indicate that qPWAWS successfully gets out of the singular point with only a small number of iterates. Besides, qPWAWS runs fast, converges with only a few iterates, obtains the same qth power median as qPWA if the latter does not get stuck, and shows a reasonable rate of linear convergence. In some cases of the online portfolio selection, the q-th power median (1 < q < 2) outperforms the 1-st power median and the 2-nd power median in both the final cumulative wealth and the Sharpe Ratio, which suggests that the q-th power median (1 < q < 2) achieves better investing performance. Therefore, the de-singularity subgradient approach is beneficial to advancing both theory and practice for the extended Weber location problem.

Future works may fall into the following aspects: 1. Find the optimal de-singularity subgradient descent in (19). 2. Establish the de-singularity subgradient theory and algorithm for the q-th power and Lp-norm median problem. 3. Extend the de-singularity subgradient theory to other related optimization problems that are also challenged by the singularity problem.

Acknowledgments

This work is supported in part by the National Natural Science Foundation of China under grants 62176103, 62276114, 62206110, in part by Guangdong Basic and Applied Basic Research Foundation under grant 2023B1515120064, in part by the Science and Technology Planning Project of Guangzhou under grants 2024A04J9896, 202206030007, Nansha District: 2023ZD001, in part by Guangdong-Macao Advanced Intelligent Computing Joint Laboratory under grant 2020B1212030003, in part by the Key Laboratory of Smart Education of Guangdong Higher Education Institutes, Jinan University under grant 2022LSYS003, and in part by the Fundamental Research Funds for the Central Universities, JNU under grant 21623202.