-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathActiveSetQP.Rmd
222 lines (186 loc) · 6.63 KB
/
ActiveSetQP.Rmd
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
---
title: "Active-Set Method for Convex QP"
output:
html_document:
toc: true
toc_depth: 4
---
# Definition of QP problem
First an object is needed that defines the QP problem with all it's constraints. Note that _I_ is just a logical vector denoting whether the constraint is an inequality (_TRUE_) or equality (_FALSE_) constraint.
```{r}
setClass('QP',
slots = c(H='matrix', c='numeric', A='matrix', b='numeric', I='logical'),
validity = function(object){TRUE})
setMethod("initialize",
'QP',
function(.Object, H, c, A, b, I, ...){
.Object <- callNextMethod()
.Object@H = H
.Object@c = c
.Object@A = A
.Object@b = b
.Object@I = I
.Object
})
QP <- function(H, c, A = matrix(), b = vector('numeric'), I = vector()){
QP <- new('QP', H = H, c = c, A = A, b = b, I = I)
QP
}
```
# Active-Set Method for Convex QP
_ActiveSet_ object as defined here is a closure. It is based on Active-Set Method for Convex QP as described in Nocedal, J. e.a., _Numerical Optimization_, 2ed, 2006, p.472. Example of how to use is further down.
Note two things here.
1. Solution of step $p_k$ is not exact, so some tolerance is needed under which components of $p_k$ are considered 0. That is the _pTol_ argument.
2. Then there is the solution of the KKT matrix. In case it is positive semi-definite, then the R base _solve()_ function can not be used. So the Moore-Penrose generalized inverse is used from _MASS_ package. _svd()_ from R base could probably also be used.
```{r}
ActiveSet <- function(QP, H = QP@H, c = QP@c, A = QP@A, b = QP@b, I = QP@I,
x0, writeLog = FALSE, pTol = 1e-10){
# Private members
a <- vector('list')
W <- vector('list')
x <- vector('list')
p <- vector('list')
kCurrent <- 0
xFinal <- NULL
# Private functions
pLog <- function(string){
if (writeLog) print(string)
}
activeConstraints <- function(xk){
# Due to numerical issues, it might be that an equality constraint
# for xk does not add up to exactly 0 in A %*% xk - b. That is
# why some tolerance is used.
acIndexes <- which(abs(A %*% xk - b) <= sqrt(.Machine$double.eps))
# All equality constraints must be in working set (always)
# Since equality constrains are never removed, it is enough
# assert that they are in the working set at the start.
if (!all(I) && !all(which(!I) %in% acIndexes)) {
stop("min$ActiveSet: Not all equality constraints active!")
}
return(acIndexes)
}
g <- function(xk){
H %*% xk + c
}
K <- function(H, A, AtMultiplier = 1){
top <- cbind(H, AtMultiplier * t(A))
bottom <- cbind(A, array(0, dim = c(nrow(A), ncol(top) - ncol(A))))
rbind(top, bottom)
}
solveKKT <- function(K, c, b){
sol <- MASS::ginv(K) %*% c(-c, b) # Moore-Penrose generalized inverse
structure(sol[1:length(c)], lambdas = sol[-(1:length(c))])
}
solvepk <- function(xk, Wk){
solveKKT(K = K(H = H, A = A[Wk,,drop=FALSE]),
c = g(xk = xk),
b = rep(0, length(Wk)))
}
pk0LagrangeMultipliers <- function(xk, Wk){
qr.solve(t(A[Wk,,drop=FALSE]), g(xk = xk))
}
pk0FoundX <- function(multipliers, Wk){
# If all multipliers >= 0 for ONLY inequality constraints,
# then that is the solution
all(multipliers[I[Wk]] >= 0)
}
pk0InequalityToRemove <- function(multipliers, Wk){
# Inequality constraint that is causing trouble.
Imultipliers <- multipliers[I[Wk]]
IWk <- Wk[I[Wk]]
IWk[which.min(Imultipliers)]
}
pk0Iteration <- function(k){
L <- pk0LagrangeMultipliers(xk = x[[k]], Wk = W[[k]])
if (pk0FoundX(multipliers = L, Wk = W[[k]])) {
pLog(paste('x* found after', kCurrent, 'steps!'))
xFinal <<- x[[k]]
return(x[[k]])
} else {
W[[k+1]] <<- setdiff(W[[k]], pk0InequalityToRemove(L, Wk = W[[k]]))
x[[k+1]] <<- x[[k]]
}
}
pkNNComputeAk <- function(xk, pk, Wk){
Imin <- intersect(setdiff(1:length(I), Wk), which(A %*% pk < 0))
alphas <- (b[Imin] - A[Imin,,drop=FALSE] %*% xk) /
(A[Imin, , drop=FALSE] %*% pk)
# We order the blocking sets according to the least alpha first.
structure(min(1, alphas),
'blockingSet' = Imin[alphas < 1][order(alphas[alphas < 1],
decreasing = FALSE)])
}
pkNNIteration <- function(k){
a[[k]] <<- pkNNComputeAk(xk = x[[k]], pk = p[[k]], Wk = W[[k]])
x[[k+1]] <<- x[[k]] + a[[k]] * as.vector(p[[k]])
if (length(attr(a[[k]], 'blockingSet') > 0)) {
# We take the blocking constrant of the lowest alpha. This is - unlike
# what Nocedal (p472) says - necessary since x_k+1 now satisfies
# the constraint of the lowest alpha and not any other. Chosing the
# wrong constraint might end up with wrong p_k+1 = 0 since the
# wrong constraint is not satisfied by x_k+1!!!
W[[k+1]] <<- union(W[[k]], attr(a[[k]], 'blockingSet')[1])
} else {
W[[k+1]] <<- W[[k]]
}
}
iterate <- function(k){
p[[k]] <<- solvepk(xk = x[[k]], Wk = W[[k]])
# Problem here is numerical precision to assume pk == 0
if (all(abs(p[[k]]) < pTol)){ # if pk == 0
pk0Iteration(k)
} else { # pk != 0
pkNNIteration(k)
}
}
nextIteration <- function(){
iterate(k = kCurrent + 1)
kCurrent <<- kCurrent + 1
}
iterateTillFound <- function(){
while (is.null(xFinal)){
nextIteration()
}
return(xFinal)
}
# Initialize
x[[1]] <- x0
W[[1]] <- activeConstraints(x[[1]])
# public functions
list(nextIteration = nextIteration,
solve = iterateTillFound,
getWorkingSetW = function(){W},
getStepP = function(){p},
getStepLengthAlpha = function(){a},
getPointsX = function(){x})
}
```
# Examples
## Example 16.2 from Nocedal (p.452-3)
```{r}
QPeEb <- QP(H = matrix(c(6,2,1,2,5,2,1,2,4), ncol = 3),
c = -c(8,3,3),
A = matrix(c(1,0,0,1,1,1), ncol = 3),
b = c(3,0),
I = c(FALSE, FALSE))
QPeEb
```
```{r}
ActiveSet(QP = QPeEb, x0 = c(3,0,0), writeLog = TRUE)$solve()
```
## Example 16.2 from Nocedal (p.475-6)
```{r}
QPeIb <- QP(H = diag(2, 2),
c = c(-2,-5),
A = matrix(data = c(+1, -2,
-1, -2,
-1, +2,
+1, +0,
+0, +1), ncol = 2, byrow = TRUE),
b = c(-2, -6, -2, 0, 0),
I = c(TRUE, TRUE, TRUE, TRUE, TRUE))
QPeIb
```
```{r}
ActiveSet(QP = QPeIb, x0 = c(2,0), writeLog = TRUE)$solve()
```