|
33 | 33 | return a + delta; |
34 | 34 | }; |
35 | 35 |
|
| 36 | + // need some basic operations on vectors, rather than adding a dependency, |
| 37 | + // just define here |
| 38 | + function zeros(x) { var r = new Array(x); for (var i = 0; i < x; ++i) { r[i] = 0; } return r; } |
| 39 | + function zerosM(x,y) { return zeros(x).map(function() { return zeros(y); }); } |
| 40 | + venn.zerosM = zerosM; |
| 41 | + venn.zeros = zeros; |
| 42 | + |
| 43 | + function dot(a, b) { |
| 44 | + var ret = 0; |
| 45 | + for (var i = 0; i < a.length; ++i) { |
| 46 | + ret += a[i] * b[i]; |
| 47 | + } |
| 48 | + return ret; |
| 49 | + } |
| 50 | + |
| 51 | + function norm2(a) { |
| 52 | + return Math.sqrt(dot(a, a)); |
| 53 | + } |
| 54 | + venn.norm2 = norm2; |
| 55 | + |
| 56 | + function multiplyBy(a, c) { |
| 57 | + for (var i = 0; i < a.length; ++i) { |
| 58 | + a[i] *= c; |
| 59 | + } |
| 60 | + } |
| 61 | + venn.multiplyBy = multiplyBy; |
| 62 | + |
36 | 63 | function weightedSum(ret, w1, v1, w2, v2) { |
37 | 64 | for (var j = 0; j < ret.length; ++j) { |
38 | 65 | ret[j] = w1 * v1[j] + w2 * v2[j]; |
|
161 | 188 | return {f : simplex[0].fx, |
162 | 189 | solution : simplex[0]}; |
163 | 190 | }; |
| 191 | + |
| 192 | + |
| 193 | + venn.minimizeConjugateGradient = function(f, initial, params) { |
| 194 | + // allocate all memory up front here, keep out of the loop for perfomance |
| 195 | + // reasons |
| 196 | + var current = {x: initial.slice(), fx: 0, fxprime: initial.slice()}, |
| 197 | + next = {x: initial.slice(), fx: 0, fxprime: initial.slice()}, |
| 198 | + yk = initial.slice(), |
| 199 | + pk, temp, |
| 200 | + a = 1, |
| 201 | + maxIterations; |
| 202 | + |
| 203 | + params = params || {}; |
| 204 | + maxIterations = params.maxIterations || initial.length * 5; |
| 205 | + |
| 206 | + current.fx = f(current.x, current.fxprime); |
| 207 | + pk = current.fxprime.slice(); |
| 208 | + multiplyBy(pk, -1); |
| 209 | + |
| 210 | + for (var i = 0; i < maxIterations; ++i) { |
| 211 | + if (params.history) { |
| 212 | + params.history.push({x: current.x.slice(), |
| 213 | + fx: current.fx, |
| 214 | + fxprime: current.fxprime.slice()}); |
| 215 | + } |
| 216 | + |
| 217 | + a = venn.wolfeLineSearch(f, pk, current, next, a); |
| 218 | + if (!a) { |
| 219 | + // faiiled to find point that satifies wolfe conditions. |
| 220 | + // reset direction for next iteration |
| 221 | + for (var j = 0; j < pk.length; ++j) { |
| 222 | + pk[j] = -1 * current.fxprime[j]; |
| 223 | + } |
| 224 | + } else { |
| 225 | + // update direction using Polak–Ribiere CG method |
| 226 | + weightedSum(yk, 1, next.fxprime, -1, current.fxprime); |
| 227 | + |
| 228 | + var delta_k = dot(current.fxprime, current.fxprime), |
| 229 | + beta_k = Math.max(0, dot(yk, next.fxprime) / delta_k); |
| 230 | + |
| 231 | + weightedSum(pk, beta_k, pk, -1, next.fxprime); |
| 232 | + |
| 233 | + temp = current; |
| 234 | + current = next; |
| 235 | + next = temp; |
| 236 | + } |
| 237 | + |
| 238 | + if (norm2(current.fxprime) <= 1e-5) { |
| 239 | + break; |
| 240 | + } |
| 241 | + } |
| 242 | + |
| 243 | + if (params.history) { |
| 244 | + params.history.push({x: current.x.slice(), |
| 245 | + fx: current.fx, |
| 246 | + fxprime: current.fxprime.slice()}); |
| 247 | + } |
| 248 | + |
| 249 | + return current; |
| 250 | + }; |
| 251 | + |
| 252 | + var c1 = 1e-6, c2 = 0.1; |
| 253 | + |
| 254 | + /// searches along line 'pk' for a point that satifies the wolfe conditions |
| 255 | + /// See 'Numerical Optimization' by Nocedal and Wright p59-60 |
| 256 | + venn.wolfeLineSearch = function(f, pk, current, next, a) { |
| 257 | + var phi0 = current.fx, phiPrime0 = dot(current.fxprime, pk), |
| 258 | + phi = phi0, phi_old = phi0, |
| 259 | + phiPrime = phiPrime0, |
| 260 | + a0 = 0; |
| 261 | + |
| 262 | + a = a || 1; |
| 263 | + |
| 264 | + function zoom(a_lo, a_high, phi_lo) { |
| 265 | + for (var iteration = 0; iteration < 16; ++iteration) { |
| 266 | + a = (a_lo + a_high)/2; |
| 267 | + weightedSum(next.x, 1.0, current.x, a, pk); |
| 268 | + phi = next.fx = f(next.x, next.fxprime); |
| 269 | + phiPrime = dot(next.fxprime, pk); |
| 270 | + |
| 271 | + if ((phi > (phi0 + c1 * a * phiPrime0)) || |
| 272 | + (phi >= phi_lo)) { |
| 273 | + a_high = a; |
| 274 | + |
| 275 | + } else { |
| 276 | + if (Math.abs(phiPrime) <= -c2 * phiPrime0) { |
| 277 | + return a; |
| 278 | + } |
| 279 | + |
| 280 | + if (phiPrime * (a_high - a_lo) >=0) { |
| 281 | + a_high = a_lo; |
| 282 | + } |
| 283 | + |
| 284 | + a_lo = a; |
| 285 | + phi_lo = phi; |
| 286 | + } |
| 287 | + } |
| 288 | + |
| 289 | + return 0; |
| 290 | + } |
| 291 | + |
| 292 | + for (var iteration = 0; iteration < 10; ++iteration) { |
| 293 | + weightedSum(next.x, 1.0, current.x, a, pk); |
| 294 | + phi = next.fx = f(next.x, next.fxprime); |
| 295 | + phiPrime = dot(next.fxprime, pk); |
| 296 | + if ((phi > (phi0 + c1 * a * phiPrime0)) || |
| 297 | + (iteration && (phi >= phi_old))) { |
| 298 | + return zoom(a0, a, phi_old); |
| 299 | + } |
| 300 | + |
| 301 | + if (Math.abs(phiPrime) <= -c2 * phiPrime0) { |
| 302 | + return a; |
| 303 | + } |
| 304 | + |
| 305 | + if (phiPrime >= 0 ) { |
| 306 | + return zoom(a, a0, phi); |
| 307 | + } |
| 308 | + |
| 309 | + phi_old = phi; |
| 310 | + a0 = a; |
| 311 | + a *= 2; |
| 312 | + } |
| 313 | + |
| 314 | + return 0; |
| 315 | + }; |
164 | 316 | })(venn); |
0 commit comments