/
simplex.go
638 lines (590 loc) · 20.3 KB
/
simplex.go
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
// Copyright ©2016 The Gonum Authors. All rights reserved.
// Use of this source code is governed by a BSD-style
// license that can be found in the LICENSE file.
// package lp implements routines for solving linear programs.
package lp
import (
"errors"
"fmt"
"math"
"gonum.org/v1/gonum/floats"
"gonum.org/v1/gonum/mat"
)
// TODO(btracey): Could have a solver structure with an abstract factorizer. With
// this transformation the same high-level code could handle both Dense and Sparse.
// TODO(btracey): Need to improve error handling. Only want to panic if condition number inf.
// TODO(btracey): Performance enhancements. There are currently lots of linear
// solves that can be improved by doing rank-one updates. For example, the swap
// step is just a rank-one update.
// TODO(btracey): Better handling on the linear solve errors. If the condition
// number is not inf and the equation solved "well", should keep moving.
var (
ErrBland = errors.New("lp: bland: all replacements are negative or cause ill-conditioned ab")
ErrInfeasible = errors.New("lp: problem is infeasible")
ErrLinSolve = errors.New("lp: linear solve failure")
ErrUnbounded = errors.New("lp: problem is unbounded")
ErrSingular = errors.New("lp: A is singular")
ErrZeroColumn = errors.New("lp: A has a column of all zeros")
ErrZeroRow = errors.New("lp: A has a row of all zeros")
)
const badShape = "lp: size mismatch"
// TODO(btracey): Should these tolerances be part of a settings struct?
const (
// initPosTol is the tolerance on the initial condition being feasible. Strictly,
// the x should be positive, but instead it must be greater than -initPosTol.
initPosTol = 1e-13
// blandNegTol is the tolerance on the value being greater than 0 in the bland test.
blandNegTol = 1e-14
// rRoundTol is the tolerance for rounding values to zero when testing if
// constraints are met.
rRoundTol = 1e-13
// dRoundTol is the tolerance for testing if values are zero for the problem
// being unbounded.
dRoundTol = 1e-13
// phaseIZeroTol tests if the Phase I problem returned a feasible solution.
phaseIZeroTol = 1e-12
// blandZeroTol is the tolerance on testing if the bland solution can move.
blandZeroTol = 1e-12
)
// Simplex solves a linear program in standard form using Danzig's Simplex
// algorithm. The standard form of a linear program is:
// minimize cᵀ x
// s.t. A*x = b
// x >= 0 .
// The input tol sets how close to the optimal solution is found (specifically,
// when the maximal reduced cost is below tol). An error will be returned if the
// problem is infeasible or unbounded. In rare cases, numeric errors can cause
// the Simplex to fail. In this case, an error will be returned along with the
// most recently found feasible solution.
//
// The Convert function can be used to transform a general LP into standard form.
//
// The input matrix A must have at least as many columns as rows, len(c) must
// equal the number of columns of A, and len(b) must equal the number of rows of
// A or Simplex will panic. A must also have full row rank and may not contain any
// columns with all zeros, or Simplex will return an error.
//
// initialBasic can be used to set the initial set of indices for a feasible
// solution to the LP. If an initial feasible solution is not known, initialBasic
// may be nil. If initialBasic is non-nil, len(initialBasic) must equal the number
// of rows of A and must be an actual feasible solution to the LP, otherwise
// Simplex will panic.
//
// A description of the Simplex algorithm can be found in Ch. 8 of
// Strang, Gilbert. "Linear Algebra and Applications." Academic, New York (1976).
// For a detailed video introduction, see lectures 11-13 of UC Math 352
// https://www.youtube.com/watch?v=ESzYPFkY3og&index=11&list=PLh464gFUoJWOmBYla3zbZbc4nv2AXez6X.
func Simplex(c []float64, A mat.Matrix, b []float64, tol float64, initialBasic []int) (optF float64, optX []float64, err error) {
ans, x, _, err := simplex(initialBasic, c, A, b, tol)
return ans, x, err
}
func simplex(initialBasic []int, c []float64, A mat.Matrix, b []float64, tol float64) (float64, []float64, []int, error) {
err := verifyInputs(initialBasic, c, A, b)
if err != nil {
if err == ErrUnbounded {
return math.Inf(-1), nil, nil, ErrUnbounded
}
return math.NaN(), nil, nil, err
}
m, n := A.Dims()
if m == n {
// Problem is exactly constrained, perform a linear solve.
bVec := mat.NewVecDense(len(b), b)
x := make([]float64, n)
xVec := mat.NewVecDense(n, x)
err := xVec.SolveVec(A, bVec)
if err != nil {
return math.NaN(), nil, nil, ErrSingular
}
for _, v := range x {
if v < 0 {
return math.NaN(), nil, nil, ErrInfeasible
}
}
f := floats.Dot(x, c)
return f, x, nil, nil
}
// There is at least one optimal solution to the LP which is at the intersection
// to a set of constraint boundaries. For a standard form LP with m variables
// and n equality constraints, at least m-n elements of x must equal zero
// at optimality. The Simplex algorithm solves the standard-form LP by starting
// at an initial constraint vertex and successively moving to adjacent constraint
// vertices. At every vertex, the set of non-zero x values is the "basic
// feasible solution". The list of non-zero x's are maintained in basicIdxs,
// the respective columns of A are in ab, and the actual non-zero values of
// x are in xb.
//
// The LP is equality constrained such that A * x = b. This can be expanded
// to
// ab * xb + an * xn = b
// where ab are the columns of a in the basic set, and an are all of the
// other columns. Since each element of xn is zero by definition, this means
// that for all feasible solutions xb = ab^-1 * b.
//
// Before the simplex algorithm can start, an initial feasible solution must
// be found. If initialBasic is non-nil a feasible solution has been supplied.
// Otherwise the "Phase I" problem must be solved to find an initial feasible
// solution.
var basicIdxs []int // The indices of the non-zero x values.
var ab *mat.Dense // The subset of columns of A listed in basicIdxs.
var xb []float64 // The non-zero elements of x. xb = ab^-1 b
if initialBasic != nil {
// InitialBasic supplied. Panic if incorrect length or infeasible.
if len(initialBasic) != m {
panic("lp: incorrect number of initial vectors")
}
ab = mat.NewDense(m, len(initialBasic), nil)
extractColumns(ab, A, initialBasic)
xb = make([]float64, m)
err = initializeFromBasic(xb, ab, b)
if err != nil {
panic(err)
}
basicIdxs = make([]int, len(initialBasic))
copy(basicIdxs, initialBasic)
} else {
// No initial basis supplied. Solve the PhaseI problem.
basicIdxs, ab, xb, err = findInitialBasic(A, b)
if err != nil {
return math.NaN(), nil, nil, err
}
}
// basicIdxs contains the indexes for an initial feasible solution,
// ab contains the extracted columns of A, and xb contains the feasible
// solution. All x not in the basic set are 0 by construction.
// nonBasicIdx is the set of nonbasic variables.
nonBasicIdx := make([]int, 0, n-m)
inBasic := make(map[int]struct{})
for _, v := range basicIdxs {
inBasic[v] = struct{}{}
}
for i := 0; i < n; i++ {
_, ok := inBasic[i]
if !ok {
nonBasicIdx = append(nonBasicIdx, i)
}
}
// cb is the subset of c for the basic variables. an and cn
// are the equivalents to ab and cb but for the nonbasic variables.
cb := make([]float64, len(basicIdxs))
for i, idx := range basicIdxs {
cb[i] = c[idx]
}
cn := make([]float64, len(nonBasicIdx))
for i, idx := range nonBasicIdx {
cn[i] = c[idx]
}
an := mat.NewDense(m, len(nonBasicIdx), nil)
extractColumns(an, A, nonBasicIdx)
bVec := mat.NewVecDense(len(b), b)
cbVec := mat.NewVecDense(len(cb), cb)
// Temporary data needed each iteration. (Described later)
r := make([]float64, n-m)
move := make([]float64, m)
// Solve the linear program starting from the initial feasible set. This is
// the "Phase 2" problem.
//
// Algorithm:
// 1) Compute the "reduced costs" for the non-basic variables. The reduced
// costs are the lagrange multipliers of the constraints.
// r = cn - anᵀ * ab¯ᵀ * cb
// 2) If all of the reduced costs are positive, no improvement is possible,
// and the solution is optimal (xn can only increase because of
// non-negativity constraints). Otherwise, the solution can be improved and
// one element will be exchanged in the basic set.
// 3) Choose the x_n with the most negative value of r. Call this value xe.
// This variable will be swapped into the basic set.
// 4) Increase xe until the next constraint boundary is met. This will happen
// when the first element in xb becomes 0. The distance xe can increase before
// a given element in xb becomes negative can be found from
// xb = Ab^-1 b - Ab^-1 An xn
// = Ab^-1 b - Ab^-1 Ae xe
// = bhat + d x_e
// xe = bhat_i / - d_i
// where Ae is the column of A corresponding to xe.
// The constraining basic index is the first index for which this is true,
// so remove the element which is min_i (bhat_i / -d_i), assuming d_i is negative.
// If no d_i is less than 0, then the problem is unbounded.
// 5) If the new xe is 0 (that is, bhat_i == 0), then this location is at
// the intersection of several constraints. Use the Bland rule instead
// of the rule in step 4 to avoid cycling.
for {
// Compute reduced costs -- r = cn - anᵀ ab¯ᵀ cb
var tmp mat.VecDense
err = tmp.SolveVec(ab.T(), cbVec)
if err != nil {
break
}
data := make([]float64, n-m)
tmp2 := mat.NewVecDense(n-m, data)
tmp2.MulVec(an.T(), &tmp)
floats.SubTo(r, cn, data)
// Replace the most negative element in the simplex. If there are no
// negative entries then the optimal solution has been found.
minIdx := floats.MinIdx(r)
if r[minIdx] >= -tol {
break
}
for i, v := range r {
if math.Abs(v) < rRoundTol {
r[i] = 0
}
}
// Compute the moving distance.
err = computeMove(move, minIdx, A, ab, xb, nonBasicIdx)
if err != nil {
if err == ErrUnbounded {
return math.Inf(-1), nil, nil, ErrUnbounded
}
break
}
// Replace the basic index along the tightest constraint.
replace := floats.MinIdx(move)
if move[replace] <= 0 {
replace, minIdx, err = replaceBland(A, ab, xb, basicIdxs, nonBasicIdx, r, move)
if err != nil {
if err == ErrUnbounded {
return math.Inf(-1), nil, nil, ErrUnbounded
}
break
}
}
// Replace the constrained basicIdx with the newIdx.
basicIdxs[replace], nonBasicIdx[minIdx] = nonBasicIdx[minIdx], basicIdxs[replace]
cb[replace], cn[minIdx] = cn[minIdx], cb[replace]
tmpCol1 := mat.Col(nil, replace, ab)
tmpCol2 := mat.Col(nil, minIdx, an)
ab.SetCol(replace, tmpCol2)
an.SetCol(minIdx, tmpCol1)
// Compute the new xb.
xbVec := mat.NewVecDense(len(xb), xb)
err = xbVec.SolveVec(ab, bVec)
if err != nil {
break
}
}
// Found the optimum successfully or died trying. The basic variables get
// their values, and the non-basic variables are all zero.
opt := floats.Dot(cb, xb)
xopt := make([]float64, n)
for i, v := range basicIdxs {
xopt[v] = xb[i]
}
return opt, xopt, basicIdxs, err
}
// computeMove computes how far can be moved replacing each index. The results
// are stored into move.
func computeMove(move []float64, minIdx int, A mat.Matrix, ab *mat.Dense, xb []float64, nonBasicIdx []int) error {
// Find ae.
col := mat.Col(nil, nonBasicIdx[minIdx], A)
aCol := mat.NewVecDense(len(col), col)
// d = - Ab^-1 Ae
nb, _ := ab.Dims()
d := make([]float64, nb)
dVec := mat.NewVecDense(nb, d)
err := dVec.SolveVec(ab, aCol)
if err != nil {
return ErrLinSolve
}
floats.Scale(-1, d)
for i, v := range d {
if math.Abs(v) < dRoundTol {
d[i] = 0
}
}
// If no di < 0, then problem is unbounded.
if floats.Min(d) >= 0 {
return ErrUnbounded
}
// move = bhat_i / - d_i, assuming d is negative.
bHat := xb // ab^-1 b
for i, v := range d {
if v >= 0 {
move[i] = math.Inf(1)
} else {
move[i] = bHat[i] / math.Abs(v)
}
}
return nil
}
// replaceBland uses the Bland rule to find the indices to swap if the minimum
// move is 0. The indices to be swapped are replace and minIdx (following the
// nomenclature in the main routine).
func replaceBland(A mat.Matrix, ab *mat.Dense, xb []float64, basicIdxs, nonBasicIdx []int, r, move []float64) (replace, minIdx int, err error) {
m, _ := A.Dims()
// Use the traditional bland rule, except don't replace a constraint which
// causes the new ab to be singular.
for i, v := range r {
if v > -blandNegTol {
continue
}
minIdx = i
err = computeMove(move, minIdx, A, ab, xb, nonBasicIdx)
if err != nil {
// Either unbounded or something went wrong.
return -1, -1, err
}
replace = floats.MinIdx(move)
if math.Abs(move[replace]) > blandZeroTol {
// Large enough that it shouldn't be a problem
return replace, minIdx, nil
}
// Find a zero index where replacement is non-singular.
biCopy := make([]int, len(basicIdxs))
for replace, v := range move {
if v > blandZeroTol {
continue
}
copy(biCopy, basicIdxs)
biCopy[replace] = nonBasicIdx[minIdx]
abTmp := mat.NewDense(m, len(biCopy), nil)
extractColumns(abTmp, A, biCopy)
// If the condition number is reasonable, use this index.
if mat.Cond(abTmp, 1) < 1e16 {
return replace, minIdx, nil
}
}
}
return -1, -1, ErrBland
}
func verifyInputs(initialBasic []int, c []float64, A mat.Matrix, b []float64) error {
m, n := A.Dims()
if m > n {
panic("lp: more equality constraints than variables")
}
if len(c) != n {
panic("lp: c vector incorrect length")
}
if len(b) != m {
panic("lp: b vector incorrect length")
}
if len(c) != n {
panic("lp: c vector incorrect length")
}
if len(initialBasic) != 0 && len(initialBasic) != m {
panic("lp: initialBasic incorrect length")
}
// Do some sanity checks so that ab does not become singular during the
// simplex solution. If the ZeroRow checks are removed then the code for
// finding a set of linearly indepent columns must be improved.
// Check that if a row of A only has zero elements that corresponding
// element in b is zero, otherwise the problem is infeasible.
// Otherwise return ErrZeroRow.
for i := 0; i < m; i++ {
isZero := true
for j := 0; j < n; j++ {
if A.At(i, j) != 0 {
isZero = false
break
}
}
if isZero && b[i] != 0 {
// Infeasible
return ErrInfeasible
} else if isZero {
return ErrZeroRow
}
}
// Check that if a column only has zero elements that the respective C vector
// is positive (otherwise unbounded). Otherwise return ErrZeroColumn.
for j := 0; j < n; j++ {
isZero := true
for i := 0; i < m; i++ {
if A.At(i, j) != 0 {
isZero = false
break
}
}
if isZero && c[j] < 0 {
return ErrUnbounded
} else if isZero {
return ErrZeroColumn
}
}
return nil
}
// initializeFromBasic initializes the basic feasible solution given a set of
// basic indices. It extracts the columns of A specified by basicIdxs and finds
// the x values at that location. These are stored into xb.
//
// If the columns of A are not linearly independent or if the initial set is not
// feasible, an error is returned.
func initializeFromBasic(xb []float64, ab *mat.Dense, b []float64) error {
m, _ := ab.Dims()
if len(xb) != m {
panic("simplex: bad xb length")
}
xbMat := mat.NewVecDense(m, xb)
err := xbMat.SolveVec(ab, mat.NewVecDense(m, b))
if err != nil {
return errors.New("lp: subcolumns of A for supplied initial basic singular")
}
// The solve ensures that the equality constraints are met (ab * xb = b).
// Thus, the solution is feasible if and only if all of the x's are positive.
allPos := true
for _, v := range xb {
if v < -initPosTol {
allPos = false
break
}
}
if !allPos {
return errors.New("lp: supplied subcolumns not a feasible solution")
}
return nil
}
// extractColumns copies the columns specified by cols into the columns of dst.
func extractColumns(dst *mat.Dense, A mat.Matrix, cols []int) {
r, c := dst.Dims()
ra, _ := A.Dims()
if ra != r {
panic("simplex: row mismatch")
}
if c != len(cols) {
panic("simplex: column mismatch")
}
col := make([]float64, r)
for j, idx := range cols {
mat.Col(col, idx, A)
dst.SetCol(j, col)
}
}
// findInitialBasic finds an initial basic solution, and returns the basic
// indices, ab, and xb.
func findInitialBasic(A mat.Matrix, b []float64) ([]int, *mat.Dense, []float64, error) {
m, n := A.Dims()
basicIdxs := findLinearlyIndependent(A)
if len(basicIdxs) != m {
return nil, nil, nil, ErrSingular
}
// It may be that this linearly independent basis is also a feasible set. If
// so, the Phase I problem can be avoided.
ab := mat.NewDense(m, len(basicIdxs), nil)
extractColumns(ab, A, basicIdxs)
xb := make([]float64, m)
err := initializeFromBasic(xb, ab, b)
if err == nil {
return basicIdxs, ab, xb, nil
}
// This set was not feasible. Instead the "Phase I" problem must be solved
// to find an initial feasible set of basis.
//
// Method: Construct an LP whose optimal solution is a feasible solution
// to the original LP.
// 1) Introduce an artificial variable x_{n+1}.
// 2) Let x_j be the most negative element of x_b (largest constraint violation).
// 3) Add the artificial variable to A with:
// a_{n+1} = b - \sum_{i in basicIdxs} a_i + a_j
// swap j with n+1 in the basicIdxs.
// 4) Define a new LP:
// minimize x_{n+1}
// subject to [A A_{n+1}][x_1 ... x_{n+1}] = b
// x, x_{n+1} >= 0
// 5) Solve this LP. If x_{n+1} != 0, then the problem is infeasible, otherwise
// the found basis can be used as an initial basis for phase II.
//
// The extra column in Step 3 is defined such that the vector of 1s is an
// initial feasible solution.
// Find the largest constraint violator.
// Compute a_{n+1} = b - \sum{i in basicIdxs}a_i + a_j. j is in basicIDx, so
// instead just subtract the basicIdx columns that are not minIDx.
minIdx := floats.MinIdx(xb)
aX1 := make([]float64, m)
copy(aX1, b)
col := make([]float64, m)
for i, v := range basicIdxs {
if i == minIdx {
continue
}
mat.Col(col, v, A)
floats.Sub(aX1, col)
}
// Construct the new LP.
// aNew = [A, a_{n+1}]
// bNew = b
// cNew = 1 for x_{n+1}
aNew := mat.NewDense(m, n+1, nil)
aNew.Copy(A)
aNew.SetCol(n, aX1)
basicIdxs[minIdx] = n // swap minIdx with n in the basic set.
c := make([]float64, n+1)
c[n] = 1
// Solve the Phase I linear program.
_, xOpt, newBasic, err := simplex(basicIdxs, c, aNew, b, 1e-10)
if err != nil {
return nil, nil, nil, fmt.Errorf("lp: error finding feasible basis: %s", err)
}
// The original LP is infeasible if the added variable has non-zero value
// in the optimal solution to the Phase I problem.
if math.Abs(xOpt[n]) > phaseIZeroTol {
return nil, nil, nil, ErrInfeasible
}
// The basis found in Phase I is a feasible solution to the original LP if
// the added variable is not in the basis.
addedIdx := -1
for i, v := range newBasic {
if v == n {
addedIdx = i
}
xb[i] = xOpt[v]
}
if addedIdx == -1 {
extractColumns(ab, A, newBasic)
return newBasic, ab, xb, nil
}
// The value of the added variable is in the basis, but it has a zero value.
// See if exchanging another variable into the basic set finds a feasible
// solution.
basicMap := make(map[int]struct{})
for _, v := range newBasic {
basicMap[v] = struct{}{}
}
var set bool
for i := range xOpt {
if _, inBasic := basicMap[i]; inBasic {
continue
}
newBasic[addedIdx] = i
if set {
mat.Col(col, i, A)
ab.SetCol(addedIdx, col)
} else {
extractColumns(ab, A, newBasic)
set = true
}
err := initializeFromBasic(xb, ab, b)
if err == nil {
return newBasic, ab, xb, nil
}
}
return nil, nil, nil, ErrInfeasible
}
// findLinearlyIndependnt finds a set of linearly independent columns of A, and
// returns the column indexes of the linearly independent columns.
func findLinearlyIndependent(A mat.Matrix) []int {
m, n := A.Dims()
idxs := make([]int, 0, m)
columns := mat.NewDense(m, m, nil)
newCol := make([]float64, m)
// Walk in reverse order because slack variables are typically the last columns
// of A.
for i := n - 1; i >= 0; i-- {
if len(idxs) == m {
break
}
mat.Col(newCol, i, A)
columns.SetCol(len(idxs), newCol)
if len(idxs) == 0 {
// A column is linearly independent from the null set.
// If all-zero column of A are allowed, this code needs to be adjusted.
idxs = append(idxs, i)
continue
}
if mat.Cond(columns.Slice(0, m, 0, len(idxs)+1), 1) > 1e12 {
// Not linearly independent.
continue
}
idxs = append(idxs, i)
}
return idxs
}