- Email: [email protected]
A linear programming procedure based on de la Vallée Poussin's minimax estimation procedure
Computational Statistics & Data Analysis 51 (2006) 453 – 456 www.elsevier.com/locate/csda
Short Communication
A linear programming procedure based on de la Vallée Poussin’s minimax estimation procedure Richard William Farebrother∗ 11 Castle Road, Bayston Hill, Shrewsbury, UK Received 3 December 2004; received in revised form 11 October 2005; accepted 13 October 2005 Available online 2 November 2005
Abstract It is shown that de la Vallée Poussin’s 1911 procedure for the solution of linear minimax estimation problems can be adjusted to solve a class of linear programming problems. A general procedure of this type should have been accessible in the 1910s, but the historical record shows that no such procedure was developed before the work of Kantorovich, Koopmans, and Dantzig in the 1940s. © 2005 Elsevier B.V. All rights reserved. Keywords: Charles de la Vallée Poussin; Least absolute deviations; Linear programming; Minimax absolute deviation; (dual) Simplex procedure
1. Introduction In this paper, we shall show that the iterative fitting procedure employed by de la Vallée Poussin (1911) for the solution of linear minimax estimation problems can easily be adjusted to solve a class of linear optimisation problems with linear inequality constraints. We suppose that the class of problems of interest to us are of the form minimise
m
cj xj
(1)
j =1
subject to the n > m inequality constraints m
aij xj bi
for i = 1, 2, . . . , n,
(2)
j =1
where the n × m matrix A = aij has full column rank m and where, for some fixed value of h, the elements ch and a1h , a2h , . . . , anh are strictly positive. Defining the additional variable z = j cj xj , we may reformulate this problem as minimise z ∗ Tel.: +01743 873 877; fax: +44 275 4812.
E-mail address: [email protected]. 0167-9473/$ - see front matter © 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.csda.2005.10.005
(3)
454
R.W. Farebrother / Computational Statistics & Data Analysis 51 (2006) 453 – 456
subject to m
aij xj bi
for i = 1, 2, . . . , n
(4)
j =1
and m
cj xj − z = 0.
(5)
j =1 ∗ = a /c times Eq. (5) from the ith inequality (4), to obtain a second revised problem with Now, we may subtract aih ih h the same feasible region:
minimise z
(6)
subject to ∗ aih z+
aij∗ xj bi
j =h
for i = 1, 2, . . . , n
(7)
and m
cj xj − z = 0,
(8)
j =1
where for j = h, aij∗ = aij − aih cj /ch
for i = 1, 2, . . . , n.
(9)
A slight rearrangement of terms yields the third (and final) version of our problem: minimise z
(10)
subject to ⎡
z ⎣bi −
j =h
and
⎡ xh = ⎣z −
⎤ aij∗ xj ⎦
∗ aih
for i = 1, 2, . . . , n
(11)
⎤ cj xj ⎦
ch .
(12)
j =h
Temporarily ignoring Eq. (12) which only serves to define xh in terms of the other m unknowns, we find that this final version of the problem is clearly a nonsymmetric variant of the minimax absolute residual problem proposed and solved by Laplace, de Prony, Cauchy, Fourier, Chebyshev, and de la Vallée Poussin, see Farebrother (1997, 1999) for details. In principle, any one of these authors should have been able to adjust his algorithm to suit the type of problem specified here by removing the features occasioned by the specific form of the minimax absolute residual problem. But, for our purposes, the most convenient solution procedure is that offered by de la Vallée Poussin (1911). 2. Solution of the reformulated problem Arbitrarily selecting a nonsingular system of m inequality constraints from the n given in Eq. (11), say the first m inequalities, we solve these inequalities for the values of the m unknowns z and x1 , x2 , . . . , xm (excluding xh ) that
R.W. Farebrother / Computational Statistics & Data Analysis 51 (2006) 453 – 456
455
satisfy the selected system of equations as equalities. If required, the corresponding value of xh may then be obtained from an application of Eq. (12). Having obtained values of x1 , x2 , . . . , xm (possibly excluding xh ) in this way, we substitute these values into the n inequalities (11) to establish whether any one of the n − m unused constraints requires a larger value of z. If so, then we have to augment the current defining set by one of the equations associated with a larger value of z (presumably either the first or the largest such value of z) and solve a similar system of m equations in m unknowns for each of the m new combinations of m equations implied by this augmented set. Now, at least one of these m combinations must yield a larger value of z and the algorithm continues with such a combination as the new defining set. Had de la Vallée Poussin addressed this problem, he would presumably have chosen the first set of equations to be found with this property but modern linear programmers would probably prefer to use the set corresponding to the largest such increase in the value of z. Further, since the value of z increases at each stage of the procedure, the associated systems of m equations are necessarily distinct, and there is only a finite number of possible combinations, so the algorithm necessarily converges to a solution of the given problem. The only difficulty arises if all the inequalities are satisfied at the first stage, for then the solution may be unbounded below and we have to check that it is not. Clearly, this variant of de la Vallée Poussin’s procedure moves from one set of m + 1 binding constraints to another. It therefore offers a simple variant of the dual simplex solution to optimisation problems of the type specified in Section 2. For further details of the proposed algorithm, see de la Vallée Poussin (1911), as outlined in Farebrother (1999, pp. 41–43) and the related work of Stiefel (1959, 1960). 3. Preliminary transformations It remains to show that linear optimisation problems of the type defined in Eqs. (1–2) can often be converted into the form required by the solution procedure of Section 3. Having deleted any zero rows as irrelevant, and having inserted superscript o’s for clarity, we may suppose that our original problem takes the form minimise
m j =1
cjo xjo
(13)
for i = 1, 2, . . . , n.
(14)
subject to m j =1
aijo xjo bi
o = co , we assume that it is possible to find a set of values w , w , . . . , w such that Setting a0j 1 2 m j m j =1
aijo wj > 0
for i = 0, 1, 2, . . . , n.
(15)
Then, normalising wh = ±1 for some suitable choice of h, we set xh = wh xho and xj = xjo − wj xh for j = h, to obtain a set of equations in the desired form, where aih = j aijo wj and aij = aijo for j = h. In passing, we note that this positivity condition cannot be satisfied if nonnegative values of y0 , y1 , y2 , . . . , yn , not all zero, exist such that n i=0
yi aijo = 0
for j = 1, 2, . . . , m.
(16)
4. Historical remarks On examining the historical record, it comes as something of a surprise to find that a solution of this type was not discovered in the early years of the twentieth century. However, potential early practitioners could have been
456
R.W. Farebrother / Computational Statistics & Data Analysis 51 (2006) 453 – 456
discouraged by the occurrence of several problems that did not satisfy the positivity condition. On the other hand, it seems reasonable to suggest that this condition would have been satisfied by a majority of the simpler problems drawn to the attention of practitioners in the early stages of their investigation. In any case, our general procedure may readily be extended to the case in which a set of values w1 , w2 , . . . , wm are found such that j aijo wj = 0 and the lower bounds in Eq. (11) are replaced by a mixture of lower and upper bounds. Indeed, Eq. (8) may be written in the form of a double inequality: m
cj xj z
j =1
m
cj x j .
(17)
j =1
Potential early practitioners may also have been discouraged by the computational complexity of the problem. But the human computers of an earlier age were not so easily dismayed by the appearance of extensive sequences of calculations. And any qualms in this regard would surely have been overcome by the prospect of financial inducements commensurate with the significant military advantage that would have accrued to any belligerent force whose strategic decisions were guided by this technique. Further, once our hypothetical early practitioners had solved a number of optimisation problems of the special type discussed here, they would surely have discovered that they did not need to perform either of the sequences of transformations outlined in Sections 2 and 4. They would thus have come close to a prototype of the familiar dual simplex procedure and would have been able to tackle problems that did not satisfy the positivity condition of Section 4. As neither of these potential problems seems to offer a sufficient reason for the late development of a practical algorithm, I am driven to the conclusion that a prototype of a simplex-like procedure was not developed in the early years of the twentieth century because problems of the type outlined in Section 2 were not of sufficient interest to civil and military planners before the 1940s. Readers will be able to form their own judgement on this question by consulting the historical papers by Schwartz (1989), Gass (1989) and Grattan-Guinness (1994).
References de la Vallée Poussin, C.J., 1911. Sur la méthode de l’approximation minimum. Ann. Soc. Sci. Bruxelles 35, 1–16 (English translation by H. E. Salger, National Bureau of Standards, USA). Farebrother, R.W., 1997. The historical development of the linear minimax absolute residual estimation procedure 1786–1960. Comput. Statist. Data Anal. 24, 455–466. Farebrother, R.W., 1999. Fitting Linear Relationships: A History of the Calculus of Observations 1750–1900. Springer, New York. Farebrother, R.W., 2002. A dual formulation of Laplace’s minimax problem. Student 4, 81–85. Gass, S., 1989. Comments on the history of linear programming. Ann. Hist. Comput. 11, 147–151. Grattan-Guinness, I., 1994. A new type of question. On the prehistory of linear and non-linear programming. In: Knoblock, E., Rowe, D. (Eds.), The History of Modern Mathematics, vol. 3. Academic Press, New York, pp. 43–89. Schwartz, B.L., 1989. The invention of linear programming. Ann. Hist. Comput. 11, 145–147. Stiefel, E., 1959. Uber diskrete und lineare Tschebyscheff-approximationen. Numer. Math. 1, 1–28. Stiefel, E., 1960. Note on Jordan elimination, linear programming and Tchebyscheff approximation. Numer. Math. 2, 1–17.
Further Reading Farebrother, R.W., 2002. A dual formulation of Laplace’s minimax problem. Student 4, 81–85.
Short Communication
A linear programming procedure based on de la Vallée Poussin’s minimax estimation procedure Richard William Farebrother∗ 11 Castle Road, Bayston Hill, Shrewsbury, UK Received 3 December 2004; received in revised form 11 October 2005; accepted 13 October 2005 Available online 2 November 2005
Abstract It is shown that de la Vallée Poussin’s 1911 procedure for the solution of linear minimax estimation problems can be adjusted to solve a class of linear programming problems. A general procedure of this type should have been accessible in the 1910s, but the historical record shows that no such procedure was developed before the work of Kantorovich, Koopmans, and Dantzig in the 1940s. © 2005 Elsevier B.V. All rights reserved. Keywords: Charles de la Vallée Poussin; Least absolute deviations; Linear programming; Minimax absolute deviation; (dual) Simplex procedure
1. Introduction In this paper, we shall show that the iterative fitting procedure employed by de la Vallée Poussin (1911) for the solution of linear minimax estimation problems can easily be adjusted to solve a class of linear optimisation problems with linear inequality constraints. We suppose that the class of problems of interest to us are of the form minimise
m
cj xj
(1)
j =1
subject to the n > m inequality constraints m
aij xj bi
for i = 1, 2, . . . , n,
(2)
j =1
where the n × m matrix A = aij has full column rank m and where, for some fixed value of h, the elements ch and a1h , a2h , . . . , anh are strictly positive. Defining the additional variable z = j cj xj , we may reformulate this problem as minimise z ∗ Tel.: +01743 873 877; fax: +44 275 4812.
E-mail address: [email protected]. 0167-9473/$ - see front matter © 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.csda.2005.10.005
(3)
454
R.W. Farebrother / Computational Statistics & Data Analysis 51 (2006) 453 – 456
subject to m
aij xj bi
for i = 1, 2, . . . , n
(4)
j =1
and m
cj xj − z = 0.
(5)
j =1 ∗ = a /c times Eq. (5) from the ith inequality (4), to obtain a second revised problem with Now, we may subtract aih ih h the same feasible region:
minimise z
(6)
subject to ∗ aih z+
aij∗ xj bi
j =h
for i = 1, 2, . . . , n
(7)
and m
cj xj − z = 0,
(8)
j =1
where for j = h, aij∗ = aij − aih cj /ch
for i = 1, 2, . . . , n.
(9)
A slight rearrangement of terms yields the third (and final) version of our problem: minimise z
(10)
subject to ⎡
z ⎣bi −
j =h
and
⎡ xh = ⎣z −
⎤ aij∗ xj ⎦
∗ aih
for i = 1, 2, . . . , n
(11)
⎤ cj xj ⎦
ch .
(12)
j =h
Temporarily ignoring Eq. (12) which only serves to define xh in terms of the other m unknowns, we find that this final version of the problem is clearly a nonsymmetric variant of the minimax absolute residual problem proposed and solved by Laplace, de Prony, Cauchy, Fourier, Chebyshev, and de la Vallée Poussin, see Farebrother (1997, 1999) for details. In principle, any one of these authors should have been able to adjust his algorithm to suit the type of problem specified here by removing the features occasioned by the specific form of the minimax absolute residual problem. But, for our purposes, the most convenient solution procedure is that offered by de la Vallée Poussin (1911). 2. Solution of the reformulated problem Arbitrarily selecting a nonsingular system of m inequality constraints from the n given in Eq. (11), say the first m inequalities, we solve these inequalities for the values of the m unknowns z and x1 , x2 , . . . , xm (excluding xh ) that
R.W. Farebrother / Computational Statistics & Data Analysis 51 (2006) 453 – 456
455
satisfy the selected system of equations as equalities. If required, the corresponding value of xh may then be obtained from an application of Eq. (12). Having obtained values of x1 , x2 , . . . , xm (possibly excluding xh ) in this way, we substitute these values into the n inequalities (11) to establish whether any one of the n − m unused constraints requires a larger value of z. If so, then we have to augment the current defining set by one of the equations associated with a larger value of z (presumably either the first or the largest such value of z) and solve a similar system of m equations in m unknowns for each of the m new combinations of m equations implied by this augmented set. Now, at least one of these m combinations must yield a larger value of z and the algorithm continues with such a combination as the new defining set. Had de la Vallée Poussin addressed this problem, he would presumably have chosen the first set of equations to be found with this property but modern linear programmers would probably prefer to use the set corresponding to the largest such increase in the value of z. Further, since the value of z increases at each stage of the procedure, the associated systems of m equations are necessarily distinct, and there is only a finite number of possible combinations, so the algorithm necessarily converges to a solution of the given problem. The only difficulty arises if all the inequalities are satisfied at the first stage, for then the solution may be unbounded below and we have to check that it is not. Clearly, this variant of de la Vallée Poussin’s procedure moves from one set of m + 1 binding constraints to another. It therefore offers a simple variant of the dual simplex solution to optimisation problems of the type specified in Section 2. For further details of the proposed algorithm, see de la Vallée Poussin (1911), as outlined in Farebrother (1999, pp. 41–43) and the related work of Stiefel (1959, 1960). 3. Preliminary transformations It remains to show that linear optimisation problems of the type defined in Eqs. (1–2) can often be converted into the form required by the solution procedure of Section 3. Having deleted any zero rows as irrelevant, and having inserted superscript o’s for clarity, we may suppose that our original problem takes the form minimise
m j =1
cjo xjo
(13)
for i = 1, 2, . . . , n.
(14)
subject to m j =1
aijo xjo bi
o = co , we assume that it is possible to find a set of values w , w , . . . , w such that Setting a0j 1 2 m j m j =1
aijo wj > 0
for i = 0, 1, 2, . . . , n.
(15)
Then, normalising wh = ±1 for some suitable choice of h, we set xh = wh xho and xj = xjo − wj xh for j = h, to obtain a set of equations in the desired form, where aih = j aijo wj and aij = aijo for j = h. In passing, we note that this positivity condition cannot be satisfied if nonnegative values of y0 , y1 , y2 , . . . , yn , not all zero, exist such that n i=0
yi aijo = 0
for j = 1, 2, . . . , m.
(16)
4. Historical remarks On examining the historical record, it comes as something of a surprise to find that a solution of this type was not discovered in the early years of the twentieth century. However, potential early practitioners could have been
456
R.W. Farebrother / Computational Statistics & Data Analysis 51 (2006) 453 – 456
discouraged by the occurrence of several problems that did not satisfy the positivity condition. On the other hand, it seems reasonable to suggest that this condition would have been satisfied by a majority of the simpler problems drawn to the attention of practitioners in the early stages of their investigation. In any case, our general procedure may readily be extended to the case in which a set of values w1 , w2 , . . . , wm are found such that j aijo wj = 0 and the lower bounds in Eq. (11) are replaced by a mixture of lower and upper bounds. Indeed, Eq. (8) may be written in the form of a double inequality: m
cj xj z
j =1
m
cj x j .
(17)
j =1
Potential early practitioners may also have been discouraged by the computational complexity of the problem. But the human computers of an earlier age were not so easily dismayed by the appearance of extensive sequences of calculations. And any qualms in this regard would surely have been overcome by the prospect of financial inducements commensurate with the significant military advantage that would have accrued to any belligerent force whose strategic decisions were guided by this technique. Further, once our hypothetical early practitioners had solved a number of optimisation problems of the special type discussed here, they would surely have discovered that they did not need to perform either of the sequences of transformations outlined in Sections 2 and 4. They would thus have come close to a prototype of the familiar dual simplex procedure and would have been able to tackle problems that did not satisfy the positivity condition of Section 4. As neither of these potential problems seems to offer a sufficient reason for the late development of a practical algorithm, I am driven to the conclusion that a prototype of a simplex-like procedure was not developed in the early years of the twentieth century because problems of the type outlined in Section 2 were not of sufficient interest to civil and military planners before the 1940s. Readers will be able to form their own judgement on this question by consulting the historical papers by Schwartz (1989), Gass (1989) and Grattan-Guinness (1994).
References de la Vallée Poussin, C.J., 1911. Sur la méthode de l’approximation minimum. Ann. Soc. Sci. Bruxelles 35, 1–16 (English translation by H. E. Salger, National Bureau of Standards, USA). Farebrother, R.W., 1997. The historical development of the linear minimax absolute residual estimation procedure 1786–1960. Comput. Statist. Data Anal. 24, 455–466. Farebrother, R.W., 1999. Fitting Linear Relationships: A History of the Calculus of Observations 1750–1900. Springer, New York. Farebrother, R.W., 2002. A dual formulation of Laplace’s minimax problem. Student 4, 81–85. Gass, S., 1989. Comments on the history of linear programming. Ann. Hist. Comput. 11, 147–151. Grattan-Guinness, I., 1994. A new type of question. On the prehistory of linear and non-linear programming. In: Knoblock, E., Rowe, D. (Eds.), The History of Modern Mathematics, vol. 3. Academic Press, New York, pp. 43–89. Schwartz, B.L., 1989. The invention of linear programming. Ann. Hist. Comput. 11, 145–147. Stiefel, E., 1959. Uber diskrete und lineare Tschebyscheff-approximationen. Numer. Math. 1, 1–28. Stiefel, E., 1960. Note on Jordan elimination, linear programming and Tchebyscheff approximation. Numer. Math. 2, 1–17.
Further Reading Farebrother, R.W., 2002. A dual formulation of Laplace’s minimax problem. Student 4, 81–85.