import math import numpy as np from numpy.linalg import inv def rosenbrock(*a): return 70 * (a[0]**2 - a[1])**2 + 5 * (a[0] - 1)**2 + 70 * (a[1]**2 - a[2])**2 + 5 * (a[1] - 1)**2 + 30 def derX1(x): return 280 * (x[0]**2 - x[1]) * x[0] + 10 * (x[1] - 1) def derX2(x): return -140 * (x[0]**2 - x[1]) + 280 * (x[1]**2 - x[2]) * x[1] + 10 * (x[1] - 1) def derX3(x): return -140 * (x[1]**2 - x[2]) def derX1X1(x): return 840 * x[0] ** 2 - 280 * x[1] def derX1X2(x): return -280 * x[0] + 10 def derX1X3(x): return 0 def derX2X2(x): return 140 + 840 * x[1] - 280 * x[2] + 10 def derX2X3(x): return -280 * x[1] def derX3X3(x): return 140 def H(x): return np.array([[derX1X1(x), derX1X2(x), derX1X3(x)], [derX1X2(x), derX2X2(x), derX2X3(x)], [derX1X3(x), derX2X3(x), derX3X3(x)] ]) def grad(x): return np.array([derX1(x), derX2(x), derX3(x)]) def norma(x): return math.sqrt(x[0]**2 + x[1]**2 + x[2]**2) def scalMul(x1, x2): return (x1 * x2).sum() def fffq(l0): a = [1., 1.] while a[-1] < l0: a.append(a[-1] + a[-2]) a.append(a[-1] + a[-2]) return a def fibchiq(a, b, epsilon, q): fib = fffq(math.fabs(b - a)/epsilon) n = len(fib) - 1 midL = a + fib[n - 2] * (b - a) / fib[n] midR = a + fib[n - 1] * (b - a) / fib[n]; while n != 2: if q(midL) <= q(midR): b = midR midR = midL midL = a + fib[n - 2] * (b - a) / fib[n] else: a = midL midL = midR midR = a + fib[n - 1] * (b - a) / fib[n] n -= 1 return (a + b)/2 def golden(a, b, e, q): A = a D = b B = A + (3 - math.sqrt(5)) * (D - A) / 2 C = A + (math.sqrt(5) - 1) * (D - A) / 2 while 1: if D - A <= 2 * e: break if q(B) <= q(C): D = C else: A = B B = A + (3 - math.sqrt(5)) * (D - A) / 2 C = A + (math.sqrt(5) - 1) * (D - A) / 2 return A + (D - A) / 2 def bibi(a, b, e, q): x1 = a x2 = b xmid = (x2 + x1) / 2. i = 2 while 1: if ((x2 - x1) <= 2 * e): break; if q(xmid - (x2 - x1) / 2. ** i) <= q(xmid + (x2 - x1) / 2. ** i) : x2 = xmid xmid = (x2 + x1) / 2. else: x1 = xmid xmid = (x2 + x1) / 2. i += 1 return (x2 + x1) / 2. #алгоритм Флетчера-Ривза def FletcherReeves(x0, e1, e2, sigma, M, func): dk_1 = -grad(x0) dk = 0 xk_1 = x0 xk = 0 xk1 = 0 wk_1 = 0 k = 0 def q(a): return rosenbrock(xk_1[0] + a * dk_1[0], xk_1[1] + a * dk_1[1], xk_1[2] + a * dk_1[2]) alphaK = func(-1000, 1000, e1 / 10000, q) xk = xk_1 + alphaK * dk_1 while 1: #print xk_1, xk, xk1 if norma(grad(xk)) <= e1: return xk, rosenbrock(xk[0], xk[1], xk[2]) if k >= M: return xk, rosenbrock(xk[0], xk[1], xk[2]) if norma(xk_1) == 0: wk_1 = 0 else: wk_1 = norma(grad(xk)) ** 2 / norma(grad(xk_1)) ** 2 dk = -grad(xk) + wk_1 * dk_1 def q(a): return rosenbrock(xk[0] + a * dk[0], xk[1] + a * dk[1], xk[2] + a * dk[2]) alphaK = func(-1000, 1000, e1 / 10000, q) xk1 = xk + alphaK * dk if norma(xk1 - xk) < sigma and math.fabs(rosenbrock(xk1[0], xk1[1], xk1[2]) - rosenbrock(xk[0], xk[1], xk[2])) < e2: return xk, rosenbrock(xk[0], xk[1], xk[2]) k += 1 xk_1 = xk xk = xk1 dk_1 = dk xxx0 = np.array([0., 0., 0.]) print "Fletcher-Reeves:" print FletcherReeves(xxx0, 0.1, 0.001, 0.001, 100, fibchiq) %timeit FletcherReeves(xxx0, 0.1, 0.001, 0.001, 25, fibchiq) #алгоритм Полака-Рибьеры def PolakRivera(x0, e1, e2, sigma, M, func, n): dk_1 = -grad(x0) dk = 0 xk_1 = x0 xk = 0 xk1 = 0 wk_1 = 0 k = 0 def q(a): return rosenbrock(xk_1[0] + a * dk_1[0], xk_1[1] + a * dk_1[1], xk_1[2] + a * dk_1[2]) alphaK = func(-1000, 1000, e1 / 10000, q) xk = xk_1 + alphaK * dk_1 while 1: #print xk_1, xk, xk1 if norma(grad(xk)) <= e1: return xk, rosenbrock(xk[0], xk[1], xk[2]) if k >= M: return xk, rosenbrock(xk[0], xk[1], xk[2]) if k % n == 0: wk_1 = 0 else: wk_1 = scalMul(grad(xk), grad(xk) - grad(xk_1))/ norma(grad(xk_1)) ** 2 #print xk_1, xk, xk1, wk_1, k dk = -grad(xk) + wk_1 * dk_1 def q(a): return rosenbrock(xk[0] + a * dk[0], xk[1] + a * dk[1], xk[2] + a * dk[2]) alphaK = func(-1000, 1000, e1 / 10000, q) xk1 = xk + alphaK * dk if norma(xk1 - xk) < sigma and math.fabs(rosenbrock(xk1[0], xk1[1], xk1[2]) - rosenbrock(xk[0], xk[1], xk[2])) < e2: return xk, rosenbrock(xk[0], xk[1], xk[2]) k += 1 xk_1 = xk xk = xk1 dk_1 = dk xxx0 = np.array([0., 0., 0.]) print "Polak-Rivera:" print PolakRivera(xxx0, 0.1, 0.001, 0.001, 25, fibchiq, 2) %timeit PolakRivera(xxx0, 0.1, 0.001, 0.001, 25, fibchiq, 2) #алгоритм Девидона-Флетчера-Пауэлла def DFP(x0, e1, e2, sigma, M, func): xk = x0 xk1 = 0 Gk = np.identity(3) Gk1 = Gk dk = -grad(xk) k = 0 while 1: if norma(grad(xk)) < e1 and k != 0: return xk, rosenbrock(xk[0], xk[1], xk[2]) if k >= M: return xk, rosenbrock(xk[0], xk[1], xk[2]) if k != 0: dgk = grad(xk1) - grad(xk) dxk = xk1 - xk if dxk.sum() == 0 or dgk.sum() == 0: s = (3, 3) dGk = np.zeros(s) else: dGk = np.dot(np.transpose(dxk), dxk) / np.dot(dxk, dgk) \ - np.dot(np.dot(np.dot(Gk, np.transpose(dgk)), dgk), np.transpose(Gk)) \ / np.dot(np.dot(dgk, Gk), np.transpose(dgk)) Gk1 = Gk + dGk dk = - np.dot(Gk1, np.transpose(grad(xk))) # multiply doesn't work :() def q(a): return rosenbrock(xk[0] + a * dk[0], xk[1] + a * dk[1], xk[2] + a * dk[2]) alpha = func(-1000, 1000, e1 / 100, q) xk1 = xk + alpha * dk if norma(xk1 - xk) < sigma and math.fabs(rosenbrock(xk1[0], xk1[1], xk1[2]) - rosenbrock(xk[0], xk[1], xk[2])) < e2: return xk1, rosenbrock(xk1[0], xk1[1], xk1[2]) k += 1 xk = xk1 xxx0 = np.array([0., 0., 0.]) print "DFP:" print DFP(xxx0, 0.1, 0.001, 0.001, 80, fibchiq) %timeit DFP(xxx0, 0.1, 0.001, 0.001, 80, fibchiq) #алгоритм Левенберга-Марквардта def LM(x0, e1, e2, sigma, M): xk = x0 xk1 = 0 mk = 10 ** 4 dk = 0 k = 0 while 1: if norma(grad(xk)) < e1: return xk1, rosenbrock(xk1[0], xk1[1], xk1[2]) if k >= M: return xk1, rosenbrock(xk1[0], xk1[1], xk1[2]) dk = - np.dot(inv(H(xk) + mk * np.identity(3)), grad(xk)) xk1 = xk + dk if rosenbrock(xk1[0], xk1[1], xk1[2]) < rosenbrock(xk[0], xk[1], xk[2]): mk = mk / 2 else: mk = 2 * mk k += 1 xk = xk1 xxx0 = np.array([0., 0., 0.]) print "LM:" print LM(xxx0, 0.1, 0.001, 0.001, 80) %timeit LM(xxx0, 0.1, 0.001, 0.001, 80)