Robust Quadratic Programming for Price Optimization

1 Introduction

Product price is the most significant factor in determining sales. Better pricing strategies (i.e., pricing for multiple products) lead to higher total profits, and determination of the best price strategy is one of the most important strategic business decisions. Price optimization has been actively studied in marketing science and revenue management [Kunz and Crone, 2014] and has already achieved great success in such industries as online retail [Ferreira et al., 2015], fast fashion [Caro and Gallien, 2012], hotels [Koushik et al., 2012; Lee, 2011], and airlines [Côté et al., 2003]. Most price optimization frameworks consist of two stages. The first is demand modeling, which reveals such complicated relationships among prices and sales quantities as price elasticity of demand [Marshall, 2009], cross price elasticity (a.k.a. cannibalization) [van Ryzin and Mahajan, 1999], and the law of diminishing marginal utility [Marshall, 2009]. The second derives the best price strategy by maximizing a utility function, typically total revenue or profit. Most existing studies for multi-product price optimization employ mixed-integer programming [Caro and Gallien, 2012; Koushik et al., 2012; Lee, 2011] due to the discrete nature of individual prices.

Recently, [Ito and Fujimaki, 2017] have proposed prescriptive price optimization. They first employ recent advanced regression techniques (e.g. non-linear, sparse) to flexibly model complex relationships between demand and price with respect to multiple products, and they have shown that the profit maximization problem can be naturally expressed by means of a binary quadratic programming (BQP) problem. [Ito and Fujimaki, 2017] has proposed a relaxation method using semi-definite programming (SDP). Although the SDP method is more efficient than existing mixed-integer programming based methods, its scalability is still limited since it theoretically requires O(M 6 ) computational time, where M is the number of products. [Ito and Fujimaki, 2016] has utilized connection between submodularity and the substitutegoods property [Koushik et al., 2012] known in economics and has proposed an extremely fast proximal gradient algorithm with network-flow optimization that, under certain conditions, guarantees the global optimality of the original BQP problem.

A key drawback in previous works is its optimism with respect to its optimization. There exist two types of uncertainty in prescriptive price optimization, i.e., system noise and estimation error. System noise represents stochastic uncertainty in demand. Estimation error occurs due to the stochastic nature of machine learning, i.e., regression coefficients vary due to stochastic changes in sales records. Although previous works have not taken into account these uncertainties, it is known in the area of stochastic optimization that such uncertainties might significantly degrade the quality of "optimal" solutions due to optimistic bias in objective values [Delage and Ye, 2010; Duchi et al., 2016].

This paper proposes a novel robust quadratic method for price optimization. Our key contributions are mainly twofold. First, we prove that uncertainty in prescriptive price optimization can be represented by a matrix normal distribution when the least square estimation is employed. This provides a natural robust formulation of price optimization as conservative lower-bound maximization. Although it is common for robust optimization to take system noise into account [BenTal et al., 2009], our study is unique in the sense that we explicitly model the uncertainty that occurs in estimation made using machine learning (i.e., estimation error). Second, we propose algorithms for robust quadratic optimization consisting of sequential relaxation to a non-robust counterpart that employs a non-robust algorithm (such as [Ito and Fujimaki, 2016; 2017]) as its sub-routine. The algorithm is guaranteed to obtain a local optimum by O(log(N/δ)) calls of the nonrobust optimization oracle, which is much fewer than that of existing algorithms [Tütüncü and Koenig, 2004]. We have experimentally evaluated our robust pricing strategy for several values of conservativeness parameters, using an artificial price simulator, and are able to show its efficiency in terms of computational time and actual revenue improvement.

2 Problem Settings

2.1 Quadratic Price Optimization

Suppose we have M products and their prices and sales quantities are denoted by x = (x1, . . . , xM)⊤ ∈ X ⊆ RM and y = (y1, . . . , yM)⊤ ∈ RM , respectively, where X is a closed bounded set. Further, let us assume that the vector y of the sales quantities follows a multi-dimensional regression with non-linear basis functions as follows:

where v : X → RN is a N-dimensional non-linear transformation of x, A∗ ∈ RM×N is the true regression coefficient matrix representing cross price elasticity [Marshall, 2009] (this paper assumes N ≥ M), and ε is system noise over RM that follows a multidimensional normal distribution N(0, Σ∗ ) with zero mean and covariance matrix Σ∗ ∈ RM×M . Without loss of generality 1 , we assume that the first M elements in v(x) are linear features, i.e., v(x) = (x⊤, N -M non-linear features)⊤. The transformation v allows us to incorporate such non-linear effects as the law of diminishing marginal utility. We here omit non-price features but it is straightforward to incorporate such features as well, as pointed in [Ito and Fujimaki, 2016].

The goal of our price optimization is to maximize gross profit as defined by (x - c)⊤y, where c = (c1, . . . , cM)⊤ ∈ RM is a cost vector. Hereinafter, for notational simplicity, this paper assumes c = 0 and thus maximizes gross revenue x⊤y. The problem can then be expressed by the following quadratic programming:

2.2 Link to Prescriptive Price Optimization

This subsection shows that the prescriptive optimization problem [Ito and Fujimaki, 2016; 2017] is reduced to (2), which considers the following general sparse additive models for demand modeling, yi = ∑Mj=1 fi,j(xj) + bi, where fi,j : Pj →R is any uni-variate feature transformation function and Pi := {pi,1, pi,2, . . . , pi,L} ⊆ R is a set of price candidates. The price optimization problem can then be formulated as follows:

1 If we don't use linear feature x, we can simply fix the first M columns of A∗ to be zero.

where X := {x | xi ∈ Pi, i = 1, . . . ,M}. Again, for simplicity, we here assume c = 0, and this does not change the discussion that follows.

Without loss of generality, we can assume that there exists a set of basis vk : R → R and coefficients a′i,j,k for k = 1, 2, . . . ,K such that

for all xj ∈ Pj . The most typical and simplest setting is K = 1 and v1(xj) = xj (linear basis), but in general we can represent (4) by having every fi,j as a basis for i, j = 1, . . . ,M. Then we have y = Fv(x) + bi, where F ∈ RM×KM is defined by Fi,(k-1)M+j = a′i,j,k . for k = 1, 2, . . . ,K and n = 1, 2, . . . , N. This shows the equivalence of (3) and (2).

3 Robust Price Optimization

3.1 Robust QP under Matrix Normal Uncertainty

Since the true model A∗ cannot be obtained in practice, a natural way is to replace it by an estimator ˆ A. Suppose we have a set of training samples denoted by {xd, yd}Dd=1 . We assume the least square estimator ˆ A of A∗ defined as follows:

For n ≥ 1, let In denote the identity matrix of size n. We define matrices V and W by

Then by the polar decomposition, there exists a matrix Γ ∈ RN×D satisfying ΓΓ⊤ = IN and V = W 1/2Γ. Let Um,n denote the random matrix over Rm×n , where each entry of Um,n is independently generated by N(0, 1) and where N represents a normal distribution.

The next proposition provides the distribution of ˆ A.

Proposition 1 indicates that ˆ A follows the matrix normal distribution [Gupta and Nagar, 1999], denoted by MN(A∗ ,Σ∗ ,W), as follows:

In order to derive a formulation for robust price optimization, we first consider the confidence region of ˆ A. Let us define the estimator ˆΣ of Σ∗ and the confidence region Cλ as follows:

where U ∈ RM×N , and ‖U‖F is the Frobenius norm of U. Then we derive the robust optimization problem with matrix normal uncertainty as

where

Here λ ≥ 0 is a scalar parameter for controlling the conservativeness (a smaller λ corresponds to a more aggressive strategy).

3.2 Probabilistic Guarantees of Robustness

This section presents an asymptotic probabilistic guarantee of our robust formulation, which can be regarded as a specialization of [Duchi et al., 2016, Theorem 6]. It ensures that we can prevent optimistic bias in objective values and shows that our robust optimum value does not asymptotically over-estimate the true optimum revenue in designed probability.

Let χ2k be a random variable generated by the chi-squared distribution with degree of freedom k.

Proposition 2. Assume that the optimum solution of (2) is unique, and that the historical strategies xd are generated by some distribution. Then we have

where ˆx := argminx∈X maxQ∈C′λ v(x)⊤Qv(x).

4 Algorithms

4.1 Tight, Tractable Upper Bounds

Let us define the function f : X →R as follows:

Our goal, then, is to obtain a strategy x ∈ X which minimizes f. We can remove the maximization in f and obtain another representation as follows:

Lemma 3. It holds that

where ‖ · ‖ is the ℓ2 norm of ·.

Proof. We have

The last equality is attained when

The term ‖L1v(x)‖‖L2v(x)‖ is a product of the square roots of quartic functions w.r.t. v(x) and is not tractable for optimization. The next proposition provides a tractable upper bound of f.

Proposition 4. For any x ∈ X and γ > 0, it holds that

where

and

The equality is attained if and only if γ‖L1v(x)‖ = ‖L2v(x)‖.

Proof. For any v = v(x) ∈ RD , we have

The equality condition is γ(v⊤M1v) = (v⊤M2v)/γ, which is equivalent to γ‖L1v(x)‖ = ‖L2v(x)‖.

The following corollary shows theoretical tightness of this upper bound, i.e., the optimal solution of the original problem is achievable with an appropriately chosen γ.

Corollary 5. It holds that

On the basis of Proposition 4 and Corollary 5, the optimization problem can be re-written as follows:

4.2 Upper BoundMinimization Algorithms

Since robust optimization is essentially more difficult than its non-robust counterpart, we make the following assumption to ensure that the non-robust counterpart is solvable.

Assumption 6. For any Q′ satisfying Q′ + Q′⊤ ⪰ ˆQ + ˆQ⊤, we have a non-robust optimization oracle to solve minx∈X v(x)⊤Q′ v(x).

This assumption immediately holds for convex non-robust counterparts and also holds for certain types of non-convex ones, such as with branch-and-bound methods [Burer and Vandenbussche, 2008], minimum cut methods for unconstrained binary quadratic programming under submodularity conditions [Kolmogorov and Zabin, 2004], and prescriptive price optimization under substitute goods property conditions [Ito and Fujimaki, 2017; 2016]. Note that this assumption is required for our theoretical analysis, and we can employ an approximation algorithm in practice.

For Eq. (21) and a given ˜ γ, Q˜ γ + Q⊤˜ γ ⪰ ˆQ + ˆQ⊤ holds where Q˜ γ := ˆQ + λ(˜ γM1 + M2/˜ γ)/2, and hence the oracle in Assumption 6 is applicable. For later convenience, we formally summarize this in the following proposition 2 .

2 Even if such ˜x is not unique, the following discussion remains valid for an arbitrary choice of such ˜x, and therefore here we do not consider uniqueness.

Algorithm 1 Golden Section Search

Require: ˆQ, L1, L2, λ, α, β, δ

Initialize a = α, b = β, r = (√5 - 1)/2

while |a - b| ≥ δ do

b ← d if h(c) < h(d), and a ← c otherwise.

end while

Output ˜x := argminx∈X g(x, ˜ γ) where ˜ γ = (a + b)/2.

Algorithm 2 Coordinate Decent

Require: ˆQ, L1, L2, λ, γ0, δ

Check γ0 > 0, and initialize r ← ∞, ˜ γ ← γ0, ˜x ← arg minx g(x, ˜ γ).

while r - g(˜x, ˜ γ) > δ and ˜ γ ̸∈ {0,∞} do

r ← g(˜x, ˜ γ)

˜x ← arg minx g(x, ˜ γ)

˜ γ ← arg minγ g(˜x, γ) by (27).

end while

Output ˜x.

Proposition 7. Under Assumption 6, for a given ˜ γ ∈ (0,∞), we can obtain ˜x := argminx∈X g(x, ˜ γ) by means of a single call of the non-robust optimization oracle.

Further, on the basis of Proposition 4, for a given ˜x, the infimum of g(˜x, γ) with respect to γ is given by

The above consideration indicates that the minimization of g w.r.t. one of x or γ is tractable when the other is fixed. This paper proposes two simple algorithms to solve Eq. (26) Let us define a function h : (0,∞) → R by

For robust prescriptive price optimization, we can employ algorithms proposed in [Ito and Fujimaki, 2017; 2016] as h(γ).

Algorithm 1 is a golden section search [Kiefer, 1953] method in which we calculate h(γ) by a single call of the oracle. An advantage of Algorithm 1 is a convergence guarantee for a fixed number of calculation of h(γ) that provides theoretical validity to our upper bound minimization (see Section 4.3 for details). Algorithm 2 is a coordinate descent algorithm that alternates optimizations of g(x, γ) w.r.t. x and γ. Though this algorithm offers no guarantee on the number of iterations needed for convergence, our experiments in Section 6 show that Algorithm 2 empirically converges with only a few iterations and is, in practice, much faster than Algorithm 1. Let us emphasize that, though more advanced optimization techniques might be applicable, our empirical results in Section 6 show that these simple algorithms converge fast enough in practice.

4.3 Convergence analysis

For a range 0 < α ≤ β, let h[α,β] : [α, β] → R denote the function restricting the domain of h onto [α, β]. Further, let us define an upper bound l by l := (‖L1‖F + ‖L2‖F )maxx∈X ‖v(x)‖ ≥ maxx∈X max{‖L1v(x)‖, ‖L2v(x)‖}. The following theorem guarantees the convergence of Algorithm 1 in a fixed number of iterations.

Theorem8. Suppose that v is continuous and that there does not exist x ∈ X such that L1v(x) = L2v(x) = 0. Under Assumption 6, for any δ′ > 0, setting α = 1/β, β = 2λl2 /δ and δ = δ′ /λβ2 l2 , Algorithm 1 outputs x with the following property by O(log(λl/δ′ )) calls of the non-robust optimization oracle: there exists a local optimum x∗ ∈ X of f satisfying f (x) ≤ f (x∗ ) + δ′ .

For comparison with an existing method, as is discussed in the next subsection, we prove in the following proposition that, under a fairly natural assumption, the function f becomes convex. In this case, f has a unique minimum, and our algorithm is able to obtain the global minimum. Here S+ is a set of positive semidefinite matrices, and Q is defined by

Proposition 9. If X is a convex set, v is affine, and Q⊆ S+, then f is a convex function.

5 RelatedWork

Studies of robust optimization with ellipsoidal uncertainty were initiated independently by [Ben-Tal and Nemirovski, 1998] and [Ghaoui and Lebret, 1997]. For robust quadratic programming minx∈X x⊤Qx, [Ben-Tal et al., 2002] have considered the case in which A of Q = A⊤A has ellipsoidal uncertainty. For the case in which Q has ellipsoidal uncertainty, which is also the case this paper considers, [Halldórsson and Tütüncü, 2003] have proposed an interior point method for a general self-concordant convex-concave function given by [Nemirovski, 1999] to this robust optimization. Further detailed information on robust optimization can be found in [Ben-Tal et al., 2009; Bertsimas et al., 2011].

In contrast to [Halldórsson and Tütüncü, 2003] which is applicable only when the convexity in Proposition 9 holds, our algorithm is applicable to non-convex robust quadratic programming under the existence of a (possibly approximate) algorithm for a non-robust counterpart. If Proposition 9 holds, both Algorithm 1 and [Halldórsson and Tütüncü, 2003] are applicable. In such a case, the non-robust counterpart is convex quadratic programming, which can be solved by √N log(1/δ) Newton steps, and thus Algorithm 1 requires O(√N log(1/δ) log(N/δ)) Newton steps if ||v(x)||∞, ||L1||∞ and ||L1||∞ are bounded 3 . If N is large, this is better than the O(N log(1/δ)) Newton steps of the algorithm in [Halldórsson and Tütüncü, 2003].

6 Experiments

We employed the state-of-the-art prescriptive price optimization of [Ito and Fujimaki, 2016] to calculate h(γ) defined by Eq. (28), and compared our robust algorithm and non-robust

counterpart. Note that the results with λ = 0 are reduced to the results of the algorithm in [Ito and Fujimaki, 2016].

6.1 Evaluation on Artificial Price Simulator

Simulation Data We conducted experiments with an artificial price simulator to support our robust formulation and theoretical analysis. We applied a similar simulation data generation process as that of the original paper [Ito and Fujimaki, 2016]. The true demand (regression) model followed (1), where v(x) = (x1, x2, . . . , xM, 1)⊤ were linear features with N = M + 1. Then the true coefficient matrix A∗ was generated by a∗i,i ∼ U([-2M,-M]), a∗i,j ∼ U([0, 2]) if i ̸= j, and a∗i,N ∼ U([M/2, 3M/2]). We defined the pricing strategies as X := {0.6, 0.7, 0.8, 0.9, 1.0}M where 1.0 was the list price and 0.9 is 10%-off. The distribution of the non-diagonal entries a∗i,j for j = 1, 2, . . . ,M were chosen to balance the effect of ith price and the other prices on the sales quantity of ith product. The linear term was arranged to obtain a moderate optimum strategy of ≈ 0.8 when the nondiagonal entries were zero. The training data ({xd, yd})Dd=1 were generated by xd,i = 1.0, 0.9, 0.8, 0.7, 0.6 with probability 0.5, 0.2, 0.1, 0.1, 0.1, respectively for all i = 1, 2, . . . ,M, where the system noise in yd followed N(0, 25IM). Here xd,i has importance because it is natural to assume that the product is sold more at the list price than at discounted prices. In each experiment, we generated 10 true models randomly and 100 training datasets were generated for each dataset. The results are averages of all runs.

Effect of Robust Formulation Fig. 1 (left) demonstrates how robust solutions changed along with λ. We observed:

• The revenues were maximized with λ = 2, 3, 4 in all cases. By appropriately choosing the robustness parameter λ, we mitigate uncertainty in price optimization and achieve both better revenue and higher robustness.
• With the small training data (D = 5M), the non-robust solutions (λ = 0) result in much lower revenue than robust solutions. This is because the statistical uncertainty in estimation was large with small training data and therefore non-robust solutions easily became poor. As the size of training data increased, the difference became smaller.
• As λ increased, the standard deviation decreased (standard deviation is plotted only for D = 5M cases for visualization purposes). This indicates that the solutions became more conservative and thus more stable.
• Even with the best chosen λ value, the obtained revenue is 5-10% worse than the true maximal revenue that was obtained with the true demand model. This difference came from estimation error, which we cannot avoid regardless of the choice of optimization algorithms. This result indicates that the true optimal value is not achievable in expectation and we must always undergo a certain reduction in the expected revenue because of the estimation error.

Probabilistic Guarantee with Finite Samples We next demonstrate that the asymptotic probabilistic guarantee on the upper-bound of the objective value holds with a practical number of training samples and reasonably large λ, and show that our robust optimization mitigates the over-estimation problem. Fig. 1 (middle) shows the result for D = 100 and M = 10, 30, 50 along with λ. The vertical axis is the probability of over-estimation P(maxQ∈C′λ v(xλ)⊤Qv(xλ) ≤ v(xλ)⊤Q∗ v(xλ)), where xλ is the robust optimum strategy with λ. We observed:

• The non-robust optimization over-estimated revenue in 7080% of the cases, and the optimism of the non-robust optimization was empirically confirmed.
• Although the empirical probabilities (blue, green, red) were larger than the theoretical asymptotic value (pale blue) in this setting (D = 100), the difference rapidly decreased with increasing λ. With a reasonably large λ value (e.g. λ ≥ 3), both the absolute value of empirical probabilities and the difference between empirical and theoretical probabilities became sufficiently small in practice.

Comparison of Algorithms We compared Algorithm 1 and Algorithm 2, with a sampling algorithm of min-max optimization that is a baseline for the robust optimization algorithm. The ith step of the sampling problem solves (12) with the replacement of U by finite samples {Uj}j=1,...,i-1, and we can then add a new Ui that is a maximizer for the current optimum solution. Observe that the sampling algorithm approximates the original robust optimization problem by adding a series of quadratic constraints, and thus in this framework the efficient price optimization algorithm [Ito and Fujimaki, 2016] is not applicable. For solving this sample approximation, we used the MIQCP solver of GUROBI optimizer version 6.0, after convex approximation of ˆQ+Uj.

Fig. 1 (right) shows convergence of the objective value over algorithm iterations. The majority of computational time was devoted to the evaluation of h(γ) using non-robust optimization oracle in our two algorithms, and to the application of the MIQCP solver in the sampling algorithm. Thus we summarized this experiment along with algorithm iteration. All algorithms achieved the same objective values at convergence.

We observed that Algorithm 1 and the sampling algorithm needed almost the same number of iterations for convergence. However, while two applications of the non-robust optimization oracle in Algorithm 2 took 0.08 second per each iteration, the computational time needed for solving MIQCP grew as the number of quadratic constraints grew, and for λ = 20 it took 274 sec in the first 10 iterations. We also observed that Algorithm 2 took only a few iterations for convergence in all cases, and confirmed that Algorithm 2 was faster in practice and required only slightly more computational time than its non-robust counterpart.

6.2 Evaluation on Real Point-of-Sales Data

We applied the proposed method to real sales history of beers 4 [Ito and Fujimaki, 2017; Wang et al., 2015]. The data consisted of prices and sales quantities on 50 products over 642 days. We used the price features for the prediction, and we obtained 50 regression models by means of the least square estimator. For each product, we set lower and upper bounds on the product price as 60% and 100% of its historically maximum price. We then conducted robust optimization for λ = 0, 30, 60, 90. Table 1 shows the optimized robust pricing strategies for each λ. We observed:

4 The data has been provided by KSP-SP Co., LTD, http:// www.ksp-sp.com.

• For many products (red), non-robust prices (λ = 0) were drastically discounted from the average prices, and the prices monotonically increased optimal prices to average prices over λ. With the most conservative setting λ = 90, the optimized prices were very close to average values. These products are considered to have large price elasticity of demand. This is natural because the average prices and largely-discounted prices can be seen as, respectively, the most "tested" strategies and the most risky.
• For many different products (green), prices changed little over λ. These products are considered to have small price elasticity of demand and the list prices (average prices) are the safest and the most profitable.
• For some products (blue), non-robust prices were close to averaged prices. They once dropped to discounted prices and then returned to the averaged prices over λ. These products are mostly high-valued and therefore are considered to have large cross price elasticity. This complex

behavior of the optimized prices was caused by the price changes in the other products.

Although we were unable to evaluate these strategies in real retail stores, we have learned here that our algorithm provides a way to simulate scenarios with different risk levels, and the users can explore for the best pricing strategies on the basis of their domain expertise.

7 Summary

This paper has proposed a novel robust quadratic optimization framework for prescriptive price optimization. Our statistical observations have revealed that the uncertainty occurring in estimation in machine learning follows a matrix normal distribution, which has lead us to formulate robust quadratic programming as a conservative upper-bound minimization. Our sequential algorithms for robust quadratic programming converge fast, both practically and theoretically, and can be implemented on the basis of existing non-robust price optimization algorithms. Experimental results on both artificial and actual price data showed that our method enables users to obtain both profitable and safe price strategies in prescriptive price optimization.