Quadratic programming with transaction costs

Computers & Operations Research 35 (2008) 18 – 33 www.elsevier.com/locate/cor

Quadratic programming with transaction costs夡 Michael J. Besta,∗ , Jaroslava Hlouskovab a Professor, Department of Combinatorics and Optimization, Faculty of Mathematics, University of Waterloo, Waterloo, Ont., Canada b Research Associate, Department of Economics and Finance, Institute for Advanced Studies, Vienna, Austria

Available online 18 April 2006

Abstract We consider the problem of maximizing the mean-variance utility function of n assets. Associated with a change in an asset’s holdings from its current or target value is a transaction cost. These must be accounted for in practical problems. A straightforward way of doing so results in a 3n-dimensional optimization problem with 3n additional constraints. This higher dimensional problem is computationally expensive to solve. We present an algorithm for solving the 3n-dimensional problem by modifying an active set quadratic programming (QP) algorithm to solve the 3n-dimensional problem as an n-dimensional problem accounting for the transaction costs implicitly rather than explicitly. The method is based on deriving the optimality conditions for the higher dimensional problem solely in terms of lower dimensional quantities and requires substantially less computational effort than any active set QP algorithm applied directly on a 3n-dimensional problem. 䉷 2006 Elsevier Ltd. All rights reserved. Keywords: Quadratic programming; Active set algorithm; Portfolio optimization; Transaction costs

1. Introduction In a financial setting, holdings of assets may be bought or sold according to some criterion, typically the investor’s utility function. Here we consider the mean-variance utility function of n assets (see [123]). In response to changed data, the mean-variance utility function implies a change should be made in an investor’s holdings. In practical situations, there is a cost associated with buying or selling an asset. This is called the transaction cost. Intuitively, if a transaction cost for a particular asset is very large it may not be advantageous to change the holdings of that asset and the holdings of that asset remain at its initial target. Alternatively, if the transaction cost is quite small, it may be advantageous to make the trade and incur the transaction cost. See also [4] for an overview of transaction costs in a variety of settings. A solution for a portfolio optimization problem with linear transaction costs is given in [5]. Their model problem assumes a diagonal covariance matrix, the budget constraint and upper bounds on all assets’ holdings. We shall show that the linear transaction costs can be accounted for by solving a quadratic programming (QP) problem having 2n additional variables and 3n additional linear constraints. The contribution of this paper is to show

夡 This research was supported by the National Science and Engineering Research Council of Canada and the Institute for Advanced Studies, Vienna, Austria. ∗ Corresponding author.

0305-0548/$ - see front matter 䉷 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.cor.2006.02.013

M.J. Best, J. Hlouskova / Computers & Operations Research 35 (2008) 18 – 33

19

how to modify any active set QP algorithm to solve a 3n-dimensional problem of a special structure as an n-dimensional problem requiring substantially less computational effort. In practical problems, the number of assets, n, can be quite large: several thousand is not uncommon. The inclusion of transaction costs increases the number of problem variables from n to 3n. Being able to solve the problem solely in terms of a single n-dimensional quantities would enable the solution of large asset problems with transaction costs. This problem was addressed in [6] where a method for solving the 3n-dimensional problem by solving a series of n-dimensional optimization problems was developed. Both the algorithm in [6] and the algorithm developed in this paper account for the transaction costs implicitly rather than explicitly. The algorithm in [6] is designed for a certain class of utility functions and is based on solving a sequence of n-dimensional problems; i.e., the transaction costs are handled externally. In contrast, the algorithm presented in this paper deals only with the mean-variance (i.e., quadratic) utility function and modifies an active set QP algorithm; i.e., the transaction costs are handled internally. Active set QP methods solve certain linear equations, which are related to minimizing the objective function in the intersection of the active constraints. The coefficient matrix of this linear system is called the working matrix. For a QP problem with n variables, the size of this matrix varies from (n, n) to (2n, 2n) [7]. Updating the factors or the inverse of the working matrix requires approximately n2 to 4n2 arithmetic operations and is the dominant portion of the per iteration computational cost. Thus, we would expect an n-dimensional QP method to solve the QP problem with transaction costs in about n2 /(3n)2 = 19 of the time required by solving it directly as a 3n variable problem. In Section 2, the notation and problem formulation are introduced. In Section 3, we will revue the steps of QP AlgorithmA formulated in [7]. This is a general QP algorithm which is equivalent to many other active set QP algorithms. In Section 4, we will use the results of Section 3 to derive an algorithm for the solution of the 3n-dimensional problem solely in terms of n-dimensional quantities. Section 5 presents computational results and Section 6 gives conclusions. In the Appendix some technical results are formulated and proved. 2. Problem formulation Throughout this paper prime ( ) denotes transposition. All vectors are column vectors unless primed. The notation zi or (z)i will be used to denote the ith component of the vector z. The problem we will address is the following 3n-dimensional quadratic programming (QP) problem: minimize: c x + p  x + + q  x − + 21 x  Cx subject to: x − x + + x − = x, ˆ Ax b, x + 0, x − 0.

 (2.1)

In (2.1), p, q, xˆ and c are given n-vectors, b is a given m-vector, C is a given (n, n) positive semi-definite matrix, A = [a1 , . . . , am ] is a given (m, n) matrix and x, x + , x − are n-vectors whose optimal values are to be determined. The problem (2.1) could be solved as a 3n-dimensional problem by any QP algorithm. An important application of the type of problem (2.1) is in portfolio optimization, see [8]. In that context c = −t, p = t p˜ and q = t q, ˜ where  denotes an n-vector of expected returns on assets, C denotes their (n, n) covariance matrix and each xi represents the amount to be invested in asset i. The function t x − (1/2)x  Cx is the mean-variance utility function. xˆ denotes a target portfolio, (x + )i , (x − )i denote the amounts purchased and sold, respectively, of asset i. We assume that a target vector xˆ is given which could represent the current holdings of assets so that purchase transaction costs are measured by the amount increased from xˆ and sales transaction costs are measured by a decrease from x. ˆ Alternatively, xˆ could represent an index fund being followed such as the S&P 500. The unit purchase and sales costs for asset i are p˜ i and q˜i , respectively. The quantity p˜  x + + q˜  x − is the total transaction cost. The parameter t is a scalar parameter denoting risk tolerance. It may be desired to solve the portfolio optimization problem for a particular t (in which case it is a QP problem) or for all t 0 (in which case it is a parametric quadratic programming (PQP) problem). Finally, Ax b denote general linear constraints. In the special case of p = q 0, (2.1) can be written as   minimize: c x + ni=1 pi |xi − xˆi | + 21 x  Cx subject to: Ax b.

(2.2)

20

M.J. Best, J. Hlouskova / Computers & Operations Research 35 (2008) 18 – 33

That is, p  x + + q  x − can be written as a weighted sum of absolute values. An additional n-vector z is introduced which is required to satisfy z x − xˆ

and

z  xˆ − x.

It follows that zi |xi − xˆi |, i = 1, . . . , n. The problem (2.2) may now be solved by solving   minimize: c x + ni=1 pi zi + 21 x  Cx subject to: Ax b, z x − x, ˆ z  xˆ − x.

(2.3)

The introduction of new variables to handle absolute values is well known. See, for example, [9]. Although (2.3) does show that (2.2) can be solved as a QP problem, (2.3) has two disadvantages: the number of variables in (2.3) is twice that of (2.2) and (2.3) includes 2n additional constraints. Our purpose here is to formulate a modification of any active set QP algorithm which will solve (2.1) (and thus (2.2)) using only the n-dimensional quantity x and for any p, q satisfying p + q > 0. 3. Algorithm A: a general QP algorithm In [7] it was shown that many active set QP algorithms are equivalent in the sense that if they are started with a common point and ties are resolved in an identical manner, the algorithms would produce the identical sequence of points. This was done by formulating a general1 QP algorithm called Algorithm A. Various QP algorithms were shown to differ only in the way in which they solved the linear equations. We will derive our algorithm for (2.1) by symbolically applying Algorithm A to (2.1). We outline the steps of Algorithm A applied to min{c x +

1 2

x  Cx | Ax b},

(3.1)

where C is (n, n) symmetric positive semi-definite matrix and c is an n-vector. We can rewrite (3.1) with explicit linear constraints as min{c x + 21 x  Cx | ai x bi , i = 1, . . . , m},

(3.2)

where A = [a1 , . . . , am ], b = (b1 , . . . , bm ) . Let R = {x | Ax b} and for any x0 ∈ R, let J (x0 ) = {i | ai x0 = bi , i = 1, . . . , m}. Thus, R is the feasible region for (3.2) and J (x0 ) is the set of indices of constraints active at x0 . A key concept for Algorithm A applied to (3.2) is that of a quasi-stationary point which is defined as follows. Definition 3.1. x0 is a quasi-stationary point (QSP) for (3.2) if: (i) x0 ∈ R, (ii) x0 is an optimal solution for min{c x + 21 x  Cx | ai x = bi , all i ∈ J (x0 )}. Observe that x0 is a QSP if and only if it is feasible and there exist numbers ui , i ∈ J (x0 ) such that  −c − Cx 0 = u i ai .

(3.3)

i∈J (x0 )

Thus, x0 is a QSP if and only if x0 satisfies all of the Karush–Kuhn–Tucker (KKT) conditions for (3.2) with the possible exception of the non-negativity of the multipliers for the active inequality constraints. Algorithm A begins with an initial feasible point x0 satisfying A0 x0 =b0 where A0 is a sub-matrix of A corresponding to constraints active at x0 and b0 is the corresponding sub-vector of b. The following assumption must be satisfied. Assumption 3.1. (i) A0 has full row rank, and (ii) s  Cs > 0 for all s  = 0 and A0 s = 0. An optimal solution for (3.2) is a QSP (namely one for which the multipliers for the active inequality constraints are non-negative). Algorithm A begins with a feasible point satisfying Assumption 3.1. Then in at most n steps it constructs a QSP with a smaller objective function value. The multipliers for this are examined and if all are non-negative an optimal 1 General in the sense that it did not specify the method of solving certain linear equations.

M.J. Best, J. Hlouskova / Computers & Operations Research 35 (2008) 18 – 33

21

solution has been formed. Otherwise, the constraint with the smallest multiplier becomes inactive and Algorithm A continues to find the next QSP having a reduced objective function value. A typical iteration j of Algorithm A begins with a feasible point xj and computes a search direction sj as part of the solution of the linear equations     s −(c + Cx j ) Hj j = , (3.4) 0 vj where

 Hj =

C Aj

Aj 0

 (3.5)

and Aj is a sub-set of the rows of A corresponding to those constraints active at xj . The ordered index set Jj is used to record the indices of the active constraints associated with each row of Aj . Once sj has been determined, the maximum feasible stepsize ˆ j is computed. If ˆ j < 1 then some previously inactive constraint has become active at xj +1 = xj + ˆ j sj . Suppose the index of this constraint is l. Then Jj +1 = Jj ∪ {l}, Aj +1 = [Aj , al ] and Hj +1 is determined accordingly. The process continues by determining the next search direction from (3.4) with j replaced by j + 1. In at most n iterations, it must occur that the maximum feasible stepsize is greater than unity. Suppose for simplicity that it occurs at the present iteration j . Then xj +1 = xj + sj is a QSP for (3.2). The multipliers for the active constraints are the components of vj in the solution of (3.4). If all components of vj are non-negative, then xj +1 is optimal and the algorithm terminates. Otherwise, suppose the smallest multiplier corresponds to constraint k. Algorithm A proceeds by deleting constraint k from the active set; i.e., Jj +1 = Jj − {k} and Aj +1 is obtained from Aj by deleting the row containing the transpose of the gradient of constraint k. There are now two possibilities: either Hj +1 is non-singular or it is singular. If Hj +1 is non-singular, then Algorithm A proceeds to obtain the next search direction from (3.4) with j replaced by j + 1. If Hj +1 is singular, then sj +1 is determined as the unique search direction satisfying Aj +1 sj +1 = 0,

Cs j +1 = 0

and ak sj +1 = −1,

where ak denotes the gradient of constraint k. If ˆ j +1 , the maximum feasible stepsize for sj +1 , satisfies ˆ j +1 = +∞, then (3.2) is unbounded from below. If ˆ j +1 < + ∞, then some previously inactive constraint becomes active at xj +1 . If l is the index of this constraint, Algorithm A proceeds by setting Jj +2 = Jj +1 ∪ {l} and Aj +2 = [Aj +1 , al ]. It is shown in [7] that Hj +2 is non-singular and Algorithm A proceeds by obtaining sj +2 from (3.4) with j replaced by j + 2. The model problem (3.2) for Algorithm A contains inequality constraints only. It is very simple to extend it to account for equality constraints in the following way. Start with x0 which satisfies both the inequality and equality constraints. Put the indices of all equality constraints in J0 . In Step 3.2 of Algorithm A never drop an equality constraint from the active set. That is, in Step 3.2 set (uj +1 )k to be the smallest of the multipliers for the inequality constraints. 4. Algorithm B: a QP algorithm with transaction costs In this section we will formulate a solution algorithm for solving (2.1) using n-dimensional quantities. The major results of this section are the detailed statement of our algorithm (Algorithm B) and Theorem 4.1, which states that in a finite number of steps, Algorithm B either solves (2.1) or determines that it is unbounded from below. In order to relate quantities from a 3n-dimensional space to an n-dimensional space, it is useful to have the following definitions. Definition 4.1. Let x, y ∈ R n . (i) x + ∈ R n is the positive portion of x with respect to y if it is defined as follows. For i = 1, . . . , n  x − yi i such that xi yi , xi+ = i 0 i such that xi < yi .

(4.1)

22

M.J. Best, J. Hlouskova / Computers & Operations Research 35 (2008) 18 – 33

(ii) x − ∈ R n is the negative portion of x with respect to y if it is defined as follows. For i = 1, . . . , n  0 i such that xi > yi , xi− = yi − xi i such that xi yi .

(4.2)

Definition 4.2. If (x  , (x + ) , (x − ) ) is a feasible solution for (2.1) and satisfies the property that not both of xi+ , xi− are strictly positive for i = 1, . . . , n; i.e., xi+ xi− = 0, i = 1, . . . , n, then we call (x  , (x + ) , (x − ) ) a proper feasible solution. We now symbolically apply Algorithm A to the 3n-dimensional problem (2.1) under the assumption that p + q > 0. By treating the variables x + and x − , the constraint x − x + + x − = xˆ as well as the constraints −x + 0, −x − 0, implicitly rather than explicitly,2 we will show that all computations can be formulated in terms of n-dimensional quantities. In symbolically applying Algorithm A to (2.1), in general, non-proper iterates must be accounted for. However, the following lemma shows that if the initial point is proper then all subsequent iterates will be proper. Lemma 4.1. Let (x0 , (x0+ ) , (x0− ) ) be a proper feasible starting point for applying Algorithm A to (2.1), let Assumption 3.1 be satisfied for (2.1) and let p + q > 0. For j = 1, . . ., let (xj , (xj+ ) , (xj− ) ) be the j th iterate so obtained. Then (xj , (xj+ ) , (xj− ) ) is proper feasible solution of (2.1). Proof. Starting with proper point (x0 , (x0+ ) , (x0− ) ) , Algorithm A does not allow a constraint to become inactive until a QSP is reached. The optimality conditions for a QSP for (2.1) (the analog of (3.3) for Definition 3.1) are − c − Cx = z + A u, − p = −z − v, − q = z − w, x − x + + x − = x, ˆ Ax b, x + 0, x − 0, ui (ai x − bi ) = 0, i = 1, . . . , m, xi+ vi = 0,

xi− wi = 0,

i = 1, . . . , n.

Here v, w and z are the multipliers for the constraints x + 0, x − 0 and x − x + − x − = x, ˆ respectively. Observe that v + w = p + q. Since p + q > 0 then v + w > 0. This implies that at least one of vi and wi is strictly positive. If xi+ is strictly positive then vi = 0 and it must be that wi > 0. But then the constraint xi− will not become inactive at the next iteration. The same argument shows that if xi− > 0 then xi+ will not become inactive at the next iteration and thus (xj , (xj+ ) , (xj− ) ) is a proper feasible solution of (2.1) for all j = 1, . . . .  Suppose we begin to apply Algorithm A to (2.1) with a proper feasible solution (x0 , (x0+ ) , (x0− ) ) . At iteration j , let (xj , (xj+ ) , (xj− ) ) denote the current iterate and let (sj , (sj+ ) , (sj− ) ) be the search direction to be determined. We let Jj denote the set of indices of the constraints Ax b which are active at xj and let Aj denote the corresponding sub-matrix of rows of A. In addition, we require three index sets as follows: ⎫ Kj ⊆ {i|(xj )i = xˆi }, ⎬ {i|(xj )i  xˆi } ⊇ Kj+ ⊇ {i|(xj )i > xˆi }, (4.3) ⎭ {i|(xj )i  xˆi } ⊇ Kj− ⊇ {i|(xj )i < xˆi }. We require that Kj , Kj+ and Kj− be a partition of {1, . . . , n}; i.e., the pair-wise intersections of these index sets are null and their union is {1, . . . , n}. Taken together, the index sets Jj , Kj , Kj+ and Kj− play the role of Jj in applying Algorithm A to the n-dimensional problem (3.1). 2 We write the constraints x +  0 and x −  0 as −x +  0 and −x −  0, respectively, so that these will be “less-than or equal to” constraints as required by the model problem for Algorithm A.

M.J. Best, J. Hlouskova / Computers & Operations Research 35 (2008) 18 – 33

23

First suppose that (xj , (xj+ ) , (xj− ) ) is not a QSP. The equality constraints for the computation of sj , sj+ and sj− are Aj sj = 0,

(4.4)

−(sj+ )i = 0

for all i ∈ Kj ,

−(sj− )i = 0 for all i ∈ Kj ,

(4.5)

−(sj+ )i = 0

for all i ∈ Kj− ,

−(sj− )i = 0 for all i ∈ Kj+ ,

(4.6)

sj − sj+ + sj− = 0,

(4.7)

where the requirements of (4.6) are a consequence of requiring the active non-negativity constraints for x + and x − to remain active. Let Dj denote the coefficient matrix for the first set of constraints of (4.5). Let ei denote the ith n-dimensional unit vector. Then the columns of Dj are just −ei , all i ∈ Kj and the first set of constraints in (4.5) can be written as Dj sj+ = 0. The second set of constraints of (4.5) can be written with the same coefficient matrix as Dj sj− = 0. In a similar manner, the two sets of constraints in (4.6) can be written as Ej+ sj+ = 0 and Ej− sj− = 0, where the columns of (Ej+ ) are precisely −ei for all i ∈ Kj− and the columns of (Ej− ) are precisely −ei for all i ∈ Kj+ . The constraints (4.4)–(4.7) can now be written in matrix form as Aj sj = 0,

(4.8)

Dj sj+ = 0,

(4.9)

Dj sj− = 0,

(4.10)

Ej+ sj+ = 0,

(4.11)

Ej− sj− = 0,

(4.12)

sj − sj+ + sj− = 0.

(4.13)

Let A˜ j denote the matrix of gradients of active constraints for (2.1). Then ⎤ ⎡ 0 0 Aj ⎢ 0 Dj 0 ⎥ ⎥ ⎢ ⎢ 0 0 D j ⎥ ⎥. ⎢ ˜ Aj = ⎢ 0 E + 0 ⎥ j ⎥ ⎢ ⎣ 0 0 E− ⎦ I

(4.14)

j

−I

I

The gradient and Hessian matrix for the objective function of (2.1) at the current iterate are     c + Cx j C 0 0 and C˜ = 0 0 0 , g˜ j = p q 0 0 0

(4.15)

respectively. We now consider the equations defining the j th search direction s˜j = (sj , (sj+ ) , (sj− ) ) in applying Algorithm A to (2.1). Let u˜ j denote the vector of multipliers for the active constraints of (2.1). Then s˜j and u˜ j are the solutions of the linear equations     −g˜ j s˜j ˜ = , (4.16) H u˜ j 0 where H˜ j =

 ˜ C A˜ j

 A˜ j 0

.

(4.17)

24

M.J. Best, J. Hlouskova / Computers & Operations Research 35 (2008) 18 – 33

Partition u˜ j as u˜ j = (uj , (vj+ ) , (vj− ) , (wj+ ) , (wj− ) , zj ) ,

(4.18)

where uj , vj+ , vj− , wj+ , wj− and zj are the multipliers for the constraints (4.8)–(4.13), respectively. With this notation, the first three partitions of (4.16) yield Cs j + Aj uj + zj = −(c + Cx j ),

(4.19)

Dj vj+ + (Ej+ ) wj+ − zj = −p,

(4.20)

Dj vj− + (Ej− ) wj− + zj = −q.

(4.21)

The remaining six partitions simply restate (4.8)–(4.13). Remark 4.1. In this derivation, vj+ is a vector of dimension |Kj | with components associated with an n-vector. Their position in the n-vector is determined by means of their corresponding unit vectors in Dj . It is convenient to have an n-vector with components of vj+ placed in their natural order and zero, placed elsewhere. Let vj+ be replaced by this new vector. Although this is an abuse of notation, it makes the presentation considerably simpler. We adopt a similar convention for vj− , wj+ and wj− . With the convention of Remark 4.1, (4.20) can be expanded according to the index sets Kj , Kj+ and Kj− as follows: for all i ∈ Kj+ ,

−(zj )i = −pi

(4.22)

−(vj+ )i − (zj )i = −pi

for all i ∈ Kj ,

(4.23)

−(wj+ )i − (zj )i = −pi

for all i ∈ Kj− .

(4.24)

−(wj− )i + (zj )i = −qi

for all i ∈ Kj+ ,

(4.25)

−(vj− )i + (zj )i = −qi

for all i ∈ Kj ,

(4.26)

Similarly, (4.21) yields

(zj )i = −qi

for all i ∈ Kj− .

(4.27)

Eqs. (4.22) and (4.27) show that (zj )i is completely specified for all i ∈ Kj+ ∪ Kj− . Let zˆ j denote the vector of the remaining components of zj , namely those (zj )i for i ∈ Kj . Note that vj+ , vj− , wj+ and wj− can be determined from (4.23)–(4.26) after determining zˆ j . Vectors sj+ and sj− are determined from (4.5)–(4.7). Further, define the n-vector cˆj according to  ci + pi if i ∈ Kj+ , (4.28) cˆj ≡ cˆi (xj ) = ci if i ∈ Kj , ci − qi if i ∈ Kj− , for i = 1, . . . , n. Then from (4.8)–(4.10), (4.13), (4.19) and (4.28), sj , uj and zˆ j are the solution of the linear system     sj −(cˆj + Cx j ) , Hj uj = 0 0 zˆ j where



C Hj = Aj −Dj

Aj 0 0

−Dj 0 0

(4.29)

 .

(4.30)

M.J. Best, J. Hlouskova / Computers & Operations Research 35 (2008) 18 – 33

25

Furthermore, Hj is non-singular if the coefficient matrix H˜ j for the 3n-linear system (4.16) defined by (4.17) is nonsingular. In addition, Hj is singular if H˜ j is singular (see Lemma A.3 in the Appendix). The important conclusion here is that the key quantities sj , uj and zˆ j can be obtained by solving the much smaller linear system (4.29) rather than having to solve the larger system (4.16). Suppose now that we have solved (4.29) to obtain sj . In Step 2, Algorithm A now computes the maximum feasible stepsize ˜ j . There are three sets of inequality constraints in (2.1); Ax b, x + 0 and x − 0. The maximum feasible stepsize for each of these is    ⎧  ⎨ min bi − ai xj  all i with a  s > 0, 1 i m, i ∈ {1, . . . , m} − J , j i j ˆ j = (4.31) ai sj  ⎩  +∞ all i with ai sj 0, 1 i m, i ∈ {1, . . . , m} − Jj ,    ⎧ ⎨ min xˆi − (xj )i  all i with (s ) < 0, i ∈ K + , j i j (sj )i  ¯ j = (4.32) ⎩ + +∞ all iwith (sj )i 0, ∈ Kj ,    ⎧ ⎨ min xˆi − (xj )i  all i with (s ) > 0, i ∈ K − , j i j (sj )i  ˇ j = (4.33) ⎩ − +∞ all i with (sj )i 0, ∈ Kj , respectively. The maximum feasible stepsize for the combined constraints, ˜ j is then the minimum of ˆ j , ¯ j and ˇ j . If ˜ j < 1, then the updating in Step 3 of Algorithm A proceeds as follows. If ˜ j is determined by ˆ j in (4.31), then let l0 be the index for which the minimum occurs. Then Jj +1 = Jj ∪ {l0 } and Aj +1 = [Aj , al0 ]. The index sets Kj , Kj+ and Kj− are unchanged for the next iteration. If ˜ j is determined by ¯ j in (4.32), then let l1 be the index for which the minimum occurs. Then l1 is deleted from Kj+ to produce Kj++1 and l1 is added to Kj to give Kj +1 . In this case, Jj and Kj− , remain unchanged for the next iteration. If ˜ j is determined by ˇ j in (4.33), then let l2 be the index for which the minimum occurs. Then l2 is deleted from Kj− to produce Kj−+1 and l2 is added to Kj to give Kj +1 . In this case, Jj and Kj+ , remain unchanged for the next iteration. Note for each of the three cases for ˜ j just discussed that Hj +1 differs from Hj by the addition of a single row and column. Whatever method is used to solve the linear equations (4.29) should take advantage of this structure. We next consider the case that ˜ j 1. Algorithm A proceeds by setting j =1 and computes (xj +1 , (xj++1 ) , (xj−+1 ) ) accordingly. This new point is a QSP and Algorithm A proceeds by examining the multipliers associated with active inequality constraints. For each i ∈ Kj , the constraints (x + )i 0 and (x − )i 0 are active. Because Kj , Kj+ and Kj− is a partition of {1, . . . , n}, it follows from (4.19)–(4.21) and the definition of Dj that the multipliers for these constraints are (vj+ )i = pi − (zj )i = pi + (c + Cx j + Cs j + Aj uj )i ,

(4.34)

(vj− )i = qi + (zj )i = qi − (c + Cx j + Cs j + Aj uj )i ,

(4.35)

where i ∈ Kj and vj+ , vj− , wj+ , wj− have been modified in accordance with Remark 4.1. Constraints (4.11) require (x + )i 0 to be active for all i ∈ Kj− . Their multipliers, wj+ can be obtained from (4.24) and (4.27) (wj+ )i = pi + qi

for all i ∈ Kj− .

(4.36)

Similarly, (4.12) requires (x − )i 0 to be active for all i ∈ Kj+ . Their multipliers wj− can be obtained from (4.22) and (4.25) (wj− )i = pi + qi

for all i ∈ Kj+ .

(4.37)

Eqs. (4.36) and (4.37), together with p + q > 0 imply that wj+ 0 and wj− 0. Consequently, these constraints are never candidates to be dropped and need not be considered further. Note that this is in agreement with Lemma 4.1.

26

M.J. Best, J. Hlouskova / Computers & Operations Research 35 (2008) 18 – 33

The multipliers for the active constraints of Ax b, uj , are obtained from the solution of (4.29). Algorithm A proceeds by finding the smallest multiplier associated with an active inequality constraint for the problem. For the case at hand, the relevant multipliers are the components of uj in the solution of (4.29) and (vj+ )i and (vj− )i for all i ∈ Kj in (4.34) and (4.35), respectively. Let k0 , k1 and k2 denote the indices of the smallest elements of uj , vj+ , and vj− , respectively. Set dj = min{(uj )k0 , (vj+ )k1 , (vj− )k2 }.

(4.38)

If dj 0, then xj +1 is optimal for (2.1) and the algorithm terminates. Otherwise, the algorithm continues depending on which value the minimum occurs in (4.38). If dj =(uj )k0 , then let k denote the index of the corresponding constraint in row k0 of Aj . Then Algorithm A proceeds by deleting constraint k from the active set. Thus, Jj +1 = Jj − {k} and Aj +1 is obtained from Aj by deleting row k0 . If dj = (vj+ )k1 , then the index k1 is removed from Kj and xk1 is allowed to increase by putting k1 in Kj++1 ; i.e., Kj +1 = Kj − {k1 } and Kj++1 = Kj+ ∪ {k1 }. Dj +1 is obtained from Dj by deleting the column containing the k1 th unit vector. In this case, Jj +1 and Kj−+1 are unchanged. If dj = (vj− )k2 , then the index k2 is removed from Kj and xk2 is allowed to decrease by putting k2 in Kj−+1 ; i.e., Kj +1 = Kj − {k2 } and Kj−+1 = Kj− ∪ {k2 }. Dj +1 is obtained from Dj by deleting the column containing the k2 th unit vector. In this case, Jj +1 and Kj++1 are unchanged. H˜ j +1 is now formed from the data resulting from the last three possibilities. H˜ j +1 may or may not be singular. If it is nonsingular, Algorithm A proceeds to the next iteration and solves (4.29) with j replaced by j + 1. If H˜ j +1 is singular, then Algorithm A computes the next search direction s˜j +1 = (sj +1 , (sj++1 ) , (sj−+1 ) ) satisfying A˜ j +1 s˜j +1 = 0,

C˜ s˜j +1 = 0

and a˜  s˜j +1 = −1,

(4.39)

where a˜ is the gradient of the active constraint which was deleted from the active set of the previous iteration. From the properties of Algorithm A, s˜j +1 is unique. By definition   (ak , 0, 0) if k = k0 , a˜ = (0, −ek , 0) if k = k1 , (0, 0, −ek ) if k = k2 , where 0 is an n-vector of zeros. Expending the first two equations of (4.39) gives Aj +1 sj +1 = 0, Dj +1 sj +1 = 0, Ej++1 sj++1 = 0,

Ej−+1 sj−+1 = 0,

sj +1 − sj++1 + sj−+1 = 0, Cs j +1 = 0.

(4.40)

But then taking into account how Aj +1 was formed and the normalization of s˜j +1 in (4.39), sj satisfies Cs j +1 = 0, 

   0 Aj s = −1 , Dj j +1 0

where the “−1” is in the same position as ak in Aj , if k = k0 , or the corresponding position of −ek , if k = k1 , or k2 . But then sj +1 is the solution of ⎡ ⎤ 0   sj +1 ⎢ 0⎥ Hj uj +1 = ⎣ (4.41) ⎦, −1 zˆ j +1 0

M.J. Best, J. Hlouskova / Computers & Operations Research 35 (2008) 18 – 33

27

which necessarily has solution uj +1 = 0 and zˆ j +1 = 0. Note that the coefficient matrix for (4.41) is Hj , which is nonsingular and therefore sj +1 is uniquely determined. In giving a detailed formulation of the algorithm, it will be convenient to use the following notation, where cˆj is defined by (4.28) Fj (x) ≡ cˆj x + 21 x  Cx, gj (x) ≡ cˆj + Cx.

(4.42)

In order to initiate the algorithm, we require an x0 such that Ax 0 b, an index set J0 ⊆ {i | ai x0 = bi , 1 i m} and index sets K0 , K0+ and K0− being a partition of {1, . . . , n} and satisfying (4.3) for j = 0. Recall that we use ei to denote the ith unit vector. Let D0 be a matrix such that the columns of D0 are −ei for all i ∈ K0 and let   C A0 −D0 H0 = A0 . (4.43) 0 0 0 −D0 0 We require the following assumption be satisfied. Assumption 4.1. (i) ai , i ∈ J0 and ei , i ∈ K0 are linearly independent; i.e.,   A0 D0 has full row rank; (ii) s  Cs > 0 for all non-zero s with   A0 s = 0; D0 (iii) Eq. (2.1) has no degenerate quasi-stationary points; i.e., the gradients of those constraints active at any quasistationary point are linearly independent. Note that Assumption 4.1 implies that H0 is non-singular (see Lemma A.2 in the Appendix). Algorithm B. Model problem: problem (2.1) Begin Initialization: Let x0 be such that Assumption 4.1 is satisfied, set j = 0 and go to Step 1. Step 1: Computation of the search direction sj If Hj is nonsingular, go to Step 1.1. Otherwise go to Step 1.2. Step 1.1: Compute the Newton direction sj and the multipliers uj , zˆ j by solving the following system of linear equations     sj −gj (xj ) Hj uj = , 0 0 zˆ j where Hj is defined by (4.30) and gj (x) is defined by (4.42). Set j = 0 and go to Step 2. Step 1.2: Compute s¯j such that Aj s¯j = 0, Dj s¯j = 0 and C s¯j = 0. Set sj = −(yj s¯j )−1 s¯j , j = 1 and go to Step 2. Step 2: Computation of the stepsize j Compute ˜ j and  ˜ j = min{ˆ j , ¯ j , ˇ j } where   l such that ⎧  ⎨ bl0 − al0 xj = min bi − ai xj  all i with a  s > 0, i ∈ {1, . . . , m}, i ∈ / J j , i j ˆ j = ai sj  al0 sj ⎩ / Jj , +∞, all i with ai sj 0, i ∈ {1, . . . , m}, i ∈

28

M.J. Best, J. Hlouskova / Computers & Operations Research 35 (2008) 18 – 33

   ⎧ ⎨ xˆl1 − (xj )l1 = min xˆi − (xj )i  all i with (s ) < 0, i ∈ K + , j i j (sj )l1 (sj )i  ¯ j = ⎩ + +∞, all i with (sj )i 0, ∈ Kj ,    ⎧ ⎨ xˆl2 − (xj )l2 = min xˆi − (xj )i  all i with (s ) > 0, i ∈ K − , j i j (sj )i  (sj )l2 ˇ j = ⎩ − +∞, all i with (sj )i 0, ∈ Kj ,  l=

l0 , l1 , l2 ,

if ˜ j = ˆ j , if ˜ j = ¯ j , if ˜ j = ˇ j .

If j = 1 then if ˜ j = +∞ then (2.1) is unbounded from below and STOP else j = ˜ j and go to Step 3 endif Else j = min{˜ j , 1} and go to Step 3. Endif Step 3: Computation of xj +1 , Jj +1 , Kj +1 , Kj++1 , and Kj−+1 Compute xj +1 = xj + j sj . If j = ˜ j , go to Step 3.1. Otherwise, go to Step 3.2. Step 3.1: Set Jj +1 = Jj ∪ {l}, Kj++1 = Kj+ , Kj +1 = Kj , Kj−+1 = Kj− , if j = ˆ j , Set

Jj +1 = Jj ,

Kj++1 = Kj+ − {l},

Kj +1 = Kj ∪ {l}, Kj−+1 = Kj− ,

if j = ¯ j ,

Set

Jj +1 = Jj ,

Kj++1 = Kj+ ,

Kj +1 = Kj ∪ {l}, Kj−+1 = Kj− − {l},

if j = ˇ j ,

then compute gj +1 (xj +1 ), form Aj +1 , Ej +1 , Hj +1 , replace j with j + 1 and go to Step 1. Step 3.2: Compute the multipliers3 + (vj )i = pi + (c + Cx j + Cs j + Aj uj )i , (vj− )i = qi − (c + Cx j + Cs j + Aj uj )i , for all i ∈ Kj . From uj obtain u˜ j , the full (i.e., m-dimensional) vector of multipliers for the constraints Ax b. Compute (u˜ j )k0 = min{(u˜ j )i |1 i m}, (vj+ )k1 = min{(vj+ )i |i ∈ Kj }, (vj− )k2 = min{(vj− )i |i ∈ Kj }. Set dj = min {(u˜ j )k0 , (vj+ )k1 , (vj− )k2 }, ⎧ ⎨ k0 , k = k1 , ⎩ k2 ,

if dj = (u˜ j )k0 , if dj = (vj+ )k1 , if dj = (vj− )k2 ,

 and

yj +1 =

ak , if k = k0 , −ek , if k = k1 , ek , if k = k2 .

3 Replace v + with an n-vector whose components are those of the unmodified v + placed in the position according to their corresponding unit j j vectors in Dj and the remaining components are zero. Modify vj− in a similar way. See Remark 4.1.

M.J. Best, J. Hlouskova / Computers & Operations Research 35 (2008) 18 – 33

29

If dj +1 0 then (xj +1 , (xj++1 ) , (xj−+1 ) ) , where xj++1 and xj−+1 are positive and negative portions of xj +1 , respectively, is an optimal solution of (2.1) and STOP. Else Set Jj +1 = Jj − {k}, Kj++1 = Kj+ , Kj +1 = Kj , Kj−+1 = Kj− , if dj = (uj )k , Set

Jj +1 = Jj ,

Kj++1 = Kj+ ∪ {k},

Kj +1 = Kj − {k}, Kj−+1 = Kj− ,

if dj = (vj+ )k ,

Set

Jj +1 = Jj ,

Kj++1 = Kj+ ,

Kj +1 = Kj − {k}, Kj−+1 = Kj− ∪ {k},

if dj = (vj− )k ,

then compute gj +1 (xj +1 ), form Aj +1 , Ej +1 , Hj +1 , replace j with j + 1 and go to Step 1. Endif End Theorem 4.1. Let p + q > 0, x0 be a point satisfying Ax 0 b and let Assumption 4.1 be satisfied. Let Algorithm B be applied to (2.1) and for j = 0, 1, . . . , let xj and sj be the iterates and search directions so obtained. For j = 0, 1, . . . , let xj+ and xj− be the positive and negative portions of xj with respect to x, ˆ respectively, and let sj+ and sj− be the  positive and negative portions of sj with respect to 0, respectively. Then (xj , (xj+ ) , (xj− ) ) and (sj , (sj+ ) , (sj− ) ) are the iterates and search directions, respectively, obtained by applying Algorithm A to (2.1) provided ties are resolved in analogous ways. In addition, after a finite number of iterations Algorithm B either determines an optimal solution for (2.1) or determines that (2.1) is unbounded from below. Proof. Since Assumption 4.1 is satisfied for Algorithm B, it follows from Lemma A.2 that Assumption 3.1 is satisfied for applying Algorithm A to (2.1) with initial point (x0 , (x0+ ) , (x0− ) ) where x0+ and x0− are positive and negative portions of x0 with respect to x. ˆ From Lemma A.1, both H˜ 0 and H0 are non-singular. ˜ For j = 0, 1, . . ., let Hj be the sequence of coefficient matrices obtained by applying Algorithm A to (2.1) and let Hj be the sequence of coefficient matrices obtained by applying Algorithm B to (2.1). If H˜ j is non-singular then by Lemma A.3(i) so too is Hj and sj , uj and zˆ j are uniquely determined. By construction, (4.8)–(4.10) are satisfied. Furthermore, vj+ , vj− , wj+ and wj− are obtained from (4.22)–(4.27) following the computation of zˆ j . These quantities satisfy the linear system (4.16) with H˜ j being the coefficient matrix. Therefore, (sj , (sj+ ) , (sj− ) ) is precisely the j th search direction obtained by applying Algorithm A to (2.1). Furthermore, the computations for the maximum feasible stepsize are identical in both cases, so (xj +1 , (xj++1 ) , (xj−+1 ) ) is precisely the (j + 1)st iterate obtained by applying Algorithm A to (2.1). If H˜ j +1 is singular then by Lemma A.3(ii) so too is Hj +1 . From Algorithm A it follows that H˜ j is non-singular and furthermore from LemmaA.3(i), Hj is also non-singular. In this case,Algorithm B determines sj +1 as part of the solution of the linear system having coefficient matrix Hj and is thus uniquely determined. Furthermore, (sj +1 , (sj++1 ) , (sj−+1 ) ) satisfies all of the requirements of (4.40) and (4.41). Therefore, (sj +1 , (sj++1 ) , (sj−+1 ) ) is precisely the (j + 1)th search direction obtained by applying Algorithm A to (2.1). Furthermore, the computations for the maximum feasible stepsize are identical in both cases, so (xj +2 , (xj++2 ) , (xj−+2 ) ) is precisely the (j +2)nd iterate obtained by applying Algorithm A to (2.1). The proof of the present theorem now follows from Theorem 2 in [7].  5. Computational results In order to test Algorithm B, we programmed it to solve the following problem ⎫ minimize: −t ( x − p  x + − q  x − ) + 21 x  Cx ⎪ ⎬ subject to: x − x + + x − = x, ˆ x1 + · · · + xn = 1, d x e, ⎪ ⎭ x + 0, x − 0.

(5.1)

We compared the results with those obtained by solving the same problem with our previous algorithm developed in [6] (Algorithm BH) and with LOQO. Algorithm BH solves 3n-dimensional problem like (2.1) by solving a sequence of n-dimensional optimization problems, which also account for transaction costs implicitly rather than explicitly.

30

M.J. Best, J. Hlouskova / Computers & Operations Research 35 (2008) 18 – 33

Table 1 Execution Times for Algorithm B vs Algorithm BH vs IPM Dimension (n)

100 200 300 400 500

Execution time Algorithm B (s)

Algorithm BH (s)

IPM (s)

0.1 1.0 3.1 6.7 15.0

0.1 1.1 4.1 11.8 22.5

0.4 2.8 13.2 22.6 42.1

Algorithm BH is designed for more general objective functions than the quadratic one and transaction costs are handled externally while in Algorithm B they are handled internally. LOQO is a system developed by Vanderbei that is based on infeasible primal-dual interior point method (IPM) applied to a sequence of quadratic approximations to the given problem [10]. We use random numbers to generate data for (5.1) so that the generated problems are similar to portfolio optimization problems. We construct an (n, n) matrix Q whose elements are random numbers in the range (−2, 2) and then form C = Q Q (which is positive semi-definite and positive definite provided Q has full rank). C could then be interpreted as a covariance matrix.4 The vector of expected returns, , is composed of random numbers in the range (0.7, 1.4). We use the uniform vector with components (x) ˆ i = 1/n, i = 1, . . . , n for the target vector. The risk aversion parameter t was set to unity (t = 1). For the transaction costs, we use n-vectors p and q, all of whose components are 0.1 for p and 0.05 for q. For the lower and upper bounds, d and e, respectively, we use 0 and 0.9, respectively. The general linear constraints are taken to be only the budget constraint x1 + · · · + xn = 1. For each dimension n = 100, 200, . . . , 500, and for all three methods, we solve the problem 10 times (using different random numbers) and the figures reported in Table 1 are average execution times. LOQO reports the execution time to read the data and to perform the optimization as two separate items. Table 1 gives only the latter time. Similarly, the execution time for Algorithm B and Algorithm BH does not include the time required to construct the data. In LOQO we requested the highest possible accuracy. For both Algorithm B and Algorithm BH, an exact solution was obtained (to within computer’s accuracy). Double precision arithmetic was used throughout for these two algorithms. Algorithm B and Algorithm BH are implemented in Fortran 77 and both IPM and our methods were run on a Pentium 4, 1.8 GHz processor with 256 Mb of RAM. The execution times are summarized in Table 1, where the label “Algorithm B” refers to the algorithm formulated in this paper, the label “Algorithm BH” refers to the algorithm developed in [6] and the label “IPM” refers to the interior point method implemented in LOQO that solves 3n-dimensional problem (5.1) directly. From Table 1 we see that Algorithm B is faster than Algorithm BH, which is roughly twice as fast as the IPM for these examples. This illustrates the utility of treating transaction costs implicitly in an n-dimensional setting rather than solving the 3n-dimensional problem explicitly. 6. Conclusion We have considered the problem of maximizing the mean-variance utility function. Transaction costs must be accounted for when changing an asset’s holdings from a target value. They can be modeled by introducing an additional 2n variables giving a 3n-dimensional optimization problem. When n is large, this 3n-dimensional problem may be computationally very expensive to solve. We have developed an algorithm, which is the modification of an active set quadratic programming algorithm, to solve a 3n-dimensional problem solely in terms of n-dimensional quantities requiring substantially less computational effort. The key idea was to treat the transaction costs implicitly rather than explicitly. 4 The size of the random numbers used to generate C is not relevant as it can be accounted for in the risk aversion factor t.

M.J. Best, J. Hlouskova / Computers & Operations Research 35 (2008) 18 – 33

31

Appendix A. The following lemma is taken from [7]. It is reported here for convenience. In the context of this lemma, let A be any (k, n) matrix and let C be any (n, n) symmetric positive semi-definite matrix. Let   C A H (A) = . A 0 Lemma A.1. Suppose A has full row rank. Then H (A) is non-singular if and only if s  Cs > 0 for all non-zero s such that As = 0. Lemma A.2. Let Assumption 4.1(i), (ii) be satisfied with feasible point x0 . Then Assumption 3.1 for Algorithm A being applied to (2.1) is satisfied with (x0 , (x0+ ) , (x0− ) ) , where x0+ and x0− are positive and negative portions of x0 with respect to x, ˆ respectively. Proof. Observe that



s˜ C˜ s˜ = (s  , (s + ) , (s − ) ) 

C 0 0

0 0 0

0 0 0



s s+ s−

 = s  Cs.

(A.1)

Consequently, we need to show that s  Cs > 0 for all s  = 0 with A0 s = 0, s = s + − s − and (s + )i = 0

for all iwith (x0+ )i = 0,

(s − )i = 0

for all i with (x0− )i = 0.

But these two last equations, together with s = s + − s − imply si = 0

for all i with (x0 )i = 0,

which is just D0 s = 0. Thus, s  Cs > 0 now follows from Assumption 4.1. Consequently, we have shown    ⎫ C 0 0 s ⎪    + − + ⎪ (s , (s ) , (s ) ) 0 0 0 > 0⎪ s ⎪ ⎪ ⎪ ⎬ 0 0 0 s−  for all (s  , (s + ) , (s −) )  = 0 such that ⎪ ⎪ s ⎪ ⎪ ⎪ + ˜ ⎪ = 0. A0 s ⎭ − s

(A.2)

From Assumption 4.1, [A0 , D0 ] has full row rank. Let A˜ 0 be given by (4.14) with j = 0. We next show that A˜ 0 has full row rank.

(A.3)

Assume to the contrary that A˜ 0 does not have full row rank. Then there exist 1 , . . . , 6 , not all zero, such that  A0 1 + 6 = 0, +   (A.4) D0 2 + (E0 ) 4 − 6 = 0, D0 3 + (E0− ) 5 + 6 = 0. From the second two partitions of (A.4) D0 (2 + 3 ) + (E0+ ) 4 + (E0− ) 5 = 0. But since K0 , K0+ and K0− are a partition of {1, . . . , n}, and each of the columns of D0 , (E0+ ) and (E0+ ) is the negative of a distinct unit vector, it follows that 2 = −3 ,

4 = 0

and

5 = 0.

(A.5)

32

M.J. Best, J. Hlouskova / Computers & Operations Research 35 (2008) 18 – 33

Using this in (A.4) gives A0 1 + D0 2 = 0.

(A.6)

If 1 = 0 and 2 = 0, then from (A.5) all 1 , . . . , 6 are zero and this is then a contradiction to our assumption that A˜ 0 does not have full row rank. Therefore, there are 1 and 2 not both zero, satisfying (A.6). But this implies that [A0 , D0 ] does not have full row rank, a contradiction. The assumption that A˜ 0 does not have full row rank leads to a contradiction and is therefore false and this verifies (A.3). Eqs. (A.2), (A.3) together with Lemma A.1 imply that H˜ 0 is indeed non-singular.  Lemma A.3. Suppose Algorithm B is initiated with a point x0 for which Assumption 4.1 is satisfied and Algorithm A ˆ is initiated with (x0 , (x0+ ) , (x0− ) ) , where x0+ and x0− are the positive and negative portions of x0 with respect to x, ˜ respectively. Let Hj be defined by (4.17) for applying Algorithm A to (2.1) and Hj be defined by (4.30) for applying Algorithm B to (2.1). Then (i) if H˜ j is non-singular then Hj is non-singular and (ii) if H˜ j is singular then Hj is singular. Proof. (i) Let Aj , Dj , Ej+ , Ej− and A˜ j be as in the statement of applying Algorithm A to (2.1). Let s be non-zero such that Aj s = 0 and Dj s = 0. Let s + and s − be the positive and negative portions of s (with respect to zero; i.e., s = s + − s − ). Then by definition Ej+ s + = 0

and

Ej− s − = 0,

(A.7)

Dj s + = 0

and

Dj s − = 0.

(A.8)

From (A.7), (A.8) and the structure of A˜ j it follows that (s  , (s + ) , (s − ) ) satisfies   s A˜ j s + = 0. s− In addition, A˜ j must have full row rank, for otherwise H˜ j would be singular. Since H˜ j is non-singular, it follows from Lemma A.1 that    C 0 0 s  +  −  (s , (s ) , (s ) ) 0 0 0 s + = s  Cs > 0. s− 0 0 0 We have shown s  Cs > 0

for all non-zero s with Aj s = 0 and Dj s = 0.

(A.9)

Next suppose the rows of [Aj , Dj ] are linearly dependent. Then there exist vectors w and v, not both zero, such that Aj w = Dj v. Define the vectors 1 , 2 , 3 and 4 accordingly to 1 = w,

2 = −v,

3 = −2

and

4 = −Aj 1 .

Next observe that           0 0 I 0 Aj  0 1 + Dj 2 + 0 3 + −I 4 = 0 . Dj I 0 0 0

(A.10)

Since not both v and w can be zero, then neither all of 1 , . . . , 4 can be zero. It thus follows from (A.10) that the rows of A˜ j corresponding to its partitions 1, 2, 3 and 6, are linearly dependent. Thus, H˜ j is singular. But this is a contradiction. The assumption that the rows of [Aj , Dj ] are linearly dependent leads to a contradiction and is therefore false. We have proven that   Aj has full row rank. (A.11) Dj From (A.9), (A.11) and Lemma A.1, it now follows that Hj is nonsingular and this completes the proof of part (i).

M.J. Best, J. Hlouskova / Computers & Operations Research 35 (2008) 18 – 33

33

(ii) If H˜ j +1 is singular, then it follows from Lemma A.1 that there exist non-zero (sj +1 , (sj++1 ) , (sj−+1 ) ) such that ⎤ ⎡ ⎤   ⎡s sj +1 C 0 0 j +1 + + (A.12) 0 0 0 ⎣ sj +1 ⎦ = 0 and A˜ j +1 ⎣ sj +1 ⎦ = 0. 0 0 0 sj−+1 sj−+1 It follows from the structure of A˜ j +1 that sj +1  = 0. Then it is straightforward to show Cs j +1 = 0, so

 Hj +1

sj +1 0 0

Aj +1 sj +1 = 0

and

Dj +1 sj +1 = 0,

 =0

and thus Hj +1 is singular.

(A.13) 

References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10]

Bodie Z, Kane A, Marcus AJ. Investments. 4th ed., Irwin McGraw-Hill; 1999. Luenberger DG. Investment science. New York: Oxford University Press; 1998. Markowitz HM. Portfolio selection: efficient diversification of investment, New York: Wiley; New Haven: Yale University Press; 1959. Schattman JB. Portfolio selection under non-convex transaction costs and capital gains taxes. PhD. thesis. Rutgers Center for Operations Research, Rutgers University, USA, 2000. Best MJ, Hlouskova J. Portfolio selection and transaction costs. Computational Optimization and Applications 2003;24(1):95–116. Best MJ, Hlouskova J. An algorithm for portfolio optimization with transaction costs. Management Science 2005;51(11):1676–88. Best MJ. Equivalence of some quadratic programming algorithms. Mathematical Programming 1984;30:71–87. Best MJ, Kale JK. Quadratic programming for large-scale portfolio optimization. In: Jessica Keyes, editor. Financial services information systems. Boca Raton, FL: CRC Press; 2000. p. 513–29. Murty KG. Operations research: deterministic optimization models. Englewood Cliffs, NJ: Prentice-Hall; 1995. Vanderbei RJ. LOQO: an interior point code for quadratic programming. Optimization Methods and Software 1999;12:451–84.

Further Reading [11] [12] [13] [14] [15]

Best MJ, Ritter K. A quadratic programming algorithm. Zeitschrift fur Operations Research 1988;32(5):271–97. Dantzig GB. Linear programming and extentions. Princeton, NJ: Princeton University Press; 1963. Fletcher R. A general quadratic programming algorithm. Journal of the Institute of Mathematics and its Applications 1971;7:76–91. Mangasarian OL. Nonlinear programming. New York: McGraw-Hill; 1969. van de Panne C, Whinston A. The symmetric formulation of the simplex method for quadratic programming. Econometrica 1969;37:507–27.