Commits

Lev Panov  committed d4ba5bc

add gradientmethods.py, test python from c

  • Participants
  • Parent commits 50b625d

Comments (0)

Files changed (4)

File labs/1d_min/1d_min.cpp

 #include <math.h>
 #include <stdio.h>
 
+#include "/usr/include/python2.7/Python.h"
+
 int calls;
 
 double Function( double x )
 
 int main( void )
 {
+    Py_Initialize();
+       PyRun_SimpleString("import pylab");
+       PyRun_SimpleString("pylab.plot(range(5))");
+       PyRun_SimpleString("pylab.show()");
+       Py_Exit(0);
+       return 0;
+    
   printf("a = %6.4lf, b = %6.4lf\n", A, B);
   Test("Extremum localization", ExtremumLoc, A, B);
   //Test("dichotony", Dichotomy, A, B);

File labs/1d_min/1d_min.pro

 
 HEADERS += \
     transport.h
+
+LIBS += -lpython2.7

File labs/2d_min/gradientmethods.py

+
+from pylab import *
+
+def function(x):
+    #return exp(2.0 * x[0] * x[0] + 3.0 * x[1] * x[1]) + 1.0 * sin(0.5 * (x[0] - x[1]))
+    return 2 * x[0] * x[0] * x[0] * x[0] + x[1] * x[1] * x[1] * x[1] + 2 * x[0] * x[1] - 4
+
+
+def functionGradient(x):
+    gradient = [0.0, 0.0]
+    #gradient[0] = 4.0 * x[0] * exp(2.0 * x[0] * x[0] + 3.0 * x[1] * x[1]) + 0.5 * cos(0.5 * (x[0] - x[1]))
+    #gradient[1] = 6.0 * x[1] * exp(2.0 * x[0] * x[0] + 3.0 * x[1] * x[1]) - 0.5 * cos(0.5 * (x[0] - x[1]))
+    gradient[0] = 8 * x[0] * x[0] * x[0] + 2 * x[1]
+    gradient[1] = 4 * x[1] * x[1] * x[1] + 2 * x[0]
+    return gradient
+
+def functionHessian(x):
+    #a = 4.0 * (4.0 * x[0] * x[0] + 1.0) * exp(2.0 * x[0] * x[0] + 3.0 * x[1] * x[1]) - 0.25 * sin(0.5 * (x[0] - x[1]))
+    #b = 24.0 * x[0] * x[1] * exp(2.0 * x[0] * x[0] + 3.0 * x[1] * x[1]) + 0.25 * sin(0.5 * (x[0] - x[1]))
+    #c = 6.0 * (6.0 * x[1] * x[1] + 1.0) * exp(2.0 * x[0] * x[0] + 3.0 * x[1] * x[1]) - 0.25 * sin(0.5 * (x[0] - x[1]))
+    a = 24 * x[0] * x[0]
+    b = 2
+    c = 12 * x[1] * x[1]
+
+    hessian = [[0, 0], [0, 0]]
+    hessian[0][0] = a
+    hessian[0][1] = b
+    hessian[1][0] = b
+    hessian[1][1] = c
+    return hessian    
+
+def vectorSum(x, y):
+    return [a + b for a, b in zip(x, y)]
+
+
+def dotProduct(x, y):
+    return sum(a * b for a, b in zip(x, y))
+
+
+def vectorMultiply(x, factor):
+    return [a * factor for a in x]
+
+
+def vectorNorm(x):
+    return sqrt(dotProduct(x, x))
+
+def matrixMultiplyVector(m, x):
+    result = [0.0, 0.0]
+    result[0] = m[0][0] * x[0] + m[0][1] * x[1]
+    result[1] = m[1][0] * x[0] + m[1][1] * x[1]
+    return result    
+
+def invertMatrix(m):
+    a = m[0][0]
+    b = m[0][1]
+    c = m[1][0]
+    d = m[1][1]
+
+    det = (a * d - b * c)
+    assert det != 0
+
+    result = [[d * 1.0 / det, -b * 1.0 / det], [-c * 1.0 / det, a * 1.0 / det]]
+    return result
+
+def gradientMethodWithStepSubdivision(f, initialX, precision):
+    def _chooseStep(x, f, direction, initialStep, eps, subdivisionFactor):
+        step = initialStep
+        gradient = functionGradient(x)
+        factor = -eps * dotProduct(gradient, gradient)
+        while f(vectorSum(x, vectorMultiply(direction, step))) - f(x) > step * factor:
+            step *= subdivisionFactor
+        return step
+
+    eps = 0.5
+    subdivisionFactor = 0.5
+    initialStep = 1.0
+    points = [initialX]
+
+    iterationsNumber = 0
+    x = initialX
+    gradient = functionGradient(x)
+    antiGradient = vectorMultiply(gradient, -1.0)
+    step = _chooseStep(x, f, antiGradient, initialStep, eps, subdivisionFactor)
+
+    while abs(f(vectorSum(x, vectorMultiply(antiGradient, step))) - f(x)) > precision:
+        antiGradient = vectorMultiply(gradient, -1.0)
+        step = _chooseStep(x, f, antiGradient, initialStep, eps, subdivisionFactor)
+        x = vectorSum(x, vectorMultiply(antiGradient, step))
+        gradient = functionGradient(x)
+        points.append(x)
+        iterationsNumber += 1
+
+    return x, iterationsNumber, points
+
+
+def gradientMethodWithQuickestDescent(f, initialX, precision):
+    def _chooseStepWithGoldenSectionMethod(a, b, f, precision, currentX, currentDirection):
+        factor = 0.5 * (3.0 - sqrt(5))
+        segmentLeftBound = a
+        segmentRightBound = b
+        segmentLength = segmentRightBound - segmentLeftBound
+
+        while segmentLength >= precision:
+            y = segmentLeftBound + factor * segmentLength
+            z = segmentRightBound - factor * segmentLength
+
+            if f(vectorSum(currentX, vectorMultiply(currentDirection, y))) <= \
+               f(vectorSum(currentX, vectorMultiply(currentDirection, z))):
+                segmentRightBound = z
+            else:
+                segmentLeftBound = y
+
+            segmentLength = segmentRightBound - segmentLeftBound
+
+        return (segmentLeftBound + segmentRightBound) / 2
+
+    stepLeftBound = 0
+    stepRightBound = 5
+    points = [initialX]
+
+    iterationsNumber = 0
+    x = initialX
+    gradient = functionGradient(x)
+    antiGradient = vectorMultiply(gradient, -1.0)
+    step = _chooseStepWithGoldenSectionMethod(stepLeftBound, stepRightBound, f, precision, x, antiGradient)
+
+    while abs(f(vectorSum(x, vectorMultiply(antiGradient, step))) - f(x)) > precision:
+        antiGradient = vectorMultiply(gradient, -1.0)
+        step = _chooseStepWithGoldenSectionMethod(stepLeftBound, stepRightBound, f, precision, x, antiGradient)
+        x = vectorSum(x, vectorMultiply(antiGradient, step))
+        gradient = functionGradient(x)
+        points.append(x)
+        iterationsNumber += 1
+
+    return x, iterationsNumber, points
+
+def newtonGradientMethodWithStepSubdivision(f, initialX, precision):
+    def _chooseStep(x, f, direction, initialStep, eps, subdivisionFactor):
+        step = initialStep
+        gradient = functionGradient(x)
+        factor = eps * dotProduct(gradient, direction)
+        while f(vectorSum(x, vectorMultiply(direction, step))) - f(x) > step * factor:
+            step *= subdivisionFactor
+        return step
+
+    eps = 0.5
+    subdivisionFactor = 0.5
+    initialStep = 1.0
+    points = [initialX]
+
+    iterationsNumber = 0
+    x = initialX
+    gradient = functionGradient(x)
+    invertedHessian = invertMatrix(functionHessian(x))
+    antiGradient = vectorMultiply(gradient, -1.0)
+    direction = matrixMultiplyVector(invertedHessian, antiGradient)
+    step = _chooseStep(x, f, direction, initialStep, eps, subdivisionFactor)
+
+    while abs(f(vectorSum(x, vectorMultiply(direction, step))) - f(x)) > precision:
+        invertedHessian = invertMatrix(functionHessian(x))
+        antiGradient = vectorMultiply(gradient, -1.0)
+        direction = matrixMultiplyVector(invertedHessian, antiGradient)
+        step = _chooseStep(x, f, direction, initialStep, eps, subdivisionFactor)
+        x = vectorSum(x, vectorMultiply(direction, step))
+        gradient = functionGradient(x)
+        points.append(x)
+
+        iterationsNumber += 1
+
+    return x, iterationsNumber, points
+
+
+def proveLinearityForSubdivisionMethod(f, exactSolution, points):
+    factors = []
+
+    for i in range(len(points) - 2):
+        factor = (f(points[i + 1]) - f(exactSolution)) / (f(points[i]) - f(exactSolution))
+        factors.append(factor)
+
+    averageFactor = sum(factors) / len(factors)
+    return averageFactor, factors
+
+def proveApplicability(a, b, llr, rlr):
+    L = np.arange(llr, rlr, 0.001)
+    X = Y = np.arange(a, b, 0.1)
+    args  = [[x, y] for x, y in zip(X, Y)]
+    for l in range(0, len(L)):
+        anyStopped = False
+        for i in range(0, len(args)):
+            stop = False
+            for j in range(1, len(args)):
+                x = args[i]
+                y = args[j]
+                gradientX = functionGradient(x)
+                gradientY = functionGradient(y)
+                if vectorNorm(vectorSum(gradientX, vectorMultiply(gradientY, -1) )) > \
+                   L[l] * vectorNorm(vectorSum(x, vectorMultiply(y, -1))):
+                    stop = True
+            if stop:
+                anyStopped = True
+                break
+        if not anyStopped:
+            print 'Applicability is proved. Lipchitz constant:', L[l]
+            return True
+
+    return False
+
+def drawPlot(plotName, f, points):
+    figure(plotName, dpi = 100)
+    spectral()
+    #X = arange(-1.0, -0.4, 0.01)
+    #Y = arange(0.1, 0.85, 0.01)
+    X = arange(-1.0, 1.0, 0.01)
+    Y = arange(-1.0, 1.0, 0.01)
+    X, Y = meshgrid(X, Y)
+    Z = [f([x, y]) for x, y in zip(X, Y)]
+    contour(X, Y, Z, 75)
+
+    xOld = points[0]
+    for x in points:
+        a = [xOld[0], x[0]]
+        b = [xOld[1], x[1]]
+        plot(a, b, color='blue', linestyle='solid')
+        xOld = x
+        plot(x[0], x[1], 'b.')
+    plot(xOld[0], xOld[1], 'r.')
+
+    savefig(plotName)
+
+def solveMinimizationProblem(method, f, initialX, precision):
+    (result, iterationsNumber, points) = method(f, initialX, precision)
+    print 'Precision:', precision
+    print 'Result: (', result[0], ',', result[1], ')'
+    print 'Function value: ', function(result)
+    print 'Steps number:', iterationsNumber
+    print
+
+    return points
+
+
+if __name__ == '__main__':
+    precisions = [0.1, 0.01, 0.001]
+#   initialX = [0.3, 0.4]
+    #initialX = [-0.9, 0.5] # ok for my func
+    initialX = [-0.95, 0.95]
+
+    applicabilityIsProved = proveApplicability(-0.5, 0.5, 40, 45)
+    if not applicabilityIsProved:
+        print 'Applicability of gradient methods is not proved'
+        exit()
+    print
+
+    print 'Gradient method with step subdivision:'
+    for precision in precisions:
+        points = solveMinimizationProblem(gradientMethodWithStepSubdivision, function, initialX, precision)
+        precisionString = "%.3f" % precision
+        precisionString = precisionString.replace(".", ",")
+        drawPlot('Step_subdivision_' + precisionString, function, points)
+
+    print 'Quickest descent gradient method:'
+    for precision in precisions:
+        points = solveMinimizationProblem(gradientMethodWithQuickestDescent, function, initialX, precision)
+        precisionString = "%.3f" % precision
+        precisionString = precisionString.replace(".", ",")
+        drawPlot('Quickest_descent_' + precisionString, function, points)
+
+    print 'Newton gradient method with step subdivision:'
+    for precision in precisions:
+        points = solveMinimizationProblem(newtonGradientMethodWithStepSubdivision, function, initialX, precision)
+        precisionString = "%.3f" % precision
+        precisionString = precisionString.replace(".", ",")
+        drawPlot('Newton_' + precisionString, function, points)
+
+    print 'Rate of convergence linearity proof for subdivision method for eps = 0.001:'
+    points = gradientMethodWithStepSubdivision(function, initialX, 0.001)[2]
+    exactResult = gradientMethodWithQuickestDescent(function, initialX, 1e-10)[0]
+    (averageFactor, factors) = proveLinearityForSubdivisionMethod(function, exactResult, points)
+    print 'Factors:', factors
+    print 'Average factor:', averageFactor
+    print
+
+    #print 'Quickest descent gradient method with optimal direction:'
+    #optimalInitX = [0.3, 0.05]
+    #points = solveMinimizationProblem(gradientMethodWithQuickestDescent, function, optimalInitX, 0.001)
+    #drawPlot('Optimal direction', function, points)
+

File labs/2d_min/min_2d.cpp

 #include <string.h>
 #include <stdlib.h>
 
+//#define SASHA
 
 double eps = 1e-3;
 
     for (int i = 0; i < a; ++i)
         y *= x;
     return y;
+//    return exp(a * log(x));
+}
+
+double pow( double x, double a ) 
+{
+   return exp(a * log(x));
 }
 
 
 //    double x1_4 = x1 * x1 * x1 * x1, x2_4 = x2 * x2 * x2 * x2;
 //    return 2 * x1_4 + x2_4 + 2 * x1 * x2 - 4;
     
+#ifndef SASHA
     return 2 * pow<4>(x1) + pow<4>(x2) + 2 * x1 * x2 - 4;
+#else
+    double x[] = { x1, x2 };
+    return exp(2.0 * x[0] * x[0] + 3.0 * x[1] * x[1]) + 1.0 * sin(0.5 * (x[0] - x[1]));
+#endif
 }
 
 
 //    y[0] = 8 * x1_3 + 2 * x[1];
 //    y[1] = 4 * x2_3 + 2 * x[0];
     
+#ifndef SASHA
     y[0] = 8 * pow<3>(x[0]) + 2 * x[1];
     y[1] = 4 * pow<3>(x[1]) + 2 * x[0];
+#else    
+    y[0] = 4.0 * x[0] * exp(2.0 * x[0] * x[0] + 3.0 * x[1] * x[1]) + 0.5 * cos(0.5 * (x[0] - x[1]));
+    y[1] = 6.0 * x[1] * exp(2.0 * x[0] * x[0] + 3.0 * x[1] * x[1]) - 0.5 * cos(0.5 * (x[0] - x[1]));
+#endif
 }
 
 
 //    double a = 12 * x[0] * x[0] + 2 * x[1] * x[1],
 //           b = 4 * x[0] * x[1],
 //           c = 24 * x[1] * x[1] + 2 * x[0] * x[0];
+#ifndef SASHA
     double a = 24 * pow<2>(x[0]),
            b = 2,
            c = 12 * pow<2>(x[1]);
+#else    
+    double a = 4.0 * (4.0 * x[0] * x[0] + 1.0) * exp(2.0 * x[0] * x[0] + 3.0 * x[1] * x[1]) - 0.25 * sin(0.5 * (x[0] - x[1])),
+           b = 24.0 * x[0] * x[1] * exp(2.0 * x[0] * x[0] + 3.0 * x[1] * x[1]) + 0.25 * sin(0.5 * (x[0] - x[1])),
+           c = 6.0 * (6.0 * x[1] * x[1] + 1.0) * exp(2.0 * x[0] * x[0] + 3.0 * x[1] * x[1]) - 0.25 * sin(0.5 * (x[0] - x[1]));
+#endif
     H[0][0] = a;
     H[0][1] = H[1][0] = b;
     H[1][1] = c;
 }
 
-//double chooseStepBySubdivision( double x[], double dir[], double subdivisionFactor )
-//{
-//    double step = 1, fnew, epsilon = eps;
-//    do
-//    {
-//        step *= subdivisionFactor;
-//        fnew = Function(x[0] - step * y[0], x[1] - step * y[1]);
-////        f(vectorSum(x, vectorMultiply(direction, step)))
-//    } while (fnew - *f > -epsilon * step * norm); 
-//}
+
+double dot( double x[], double y[] )
+{
+    return x[0] * y[0] + x[1] * y[1];
+}
+
+
+void mul( double fac, double x[], double y[] )
+{
+    y[0] = fac * x[0];
+    y[1] = fac * x[1];
+}
+
+
+double chooseStepBySubdivision( double x[], double *f, double dir[], double subdivisionFactor )
+{
+    double step = 1.0, fnew, epsilon = /*eps*/0.5, grad[2], factor, xk[2];
+    
+    Gradient(x, grad);
+    factor = epsilon * dot(grad, dir);
+    
+    do
+    {
+        step *= subdivisionFactor;
+        
+        mul(step, dir, xk);
+        fnew = Function(x[0] + xk[0], x[1] + xk[1]);
+    } while (fnew - *f > step * factor); 
+    
+    return step;
+}
+
 
 void Gradient2( double x[], double y[], int recalculate )
 {
 }
 
 
+double Step_Newton( double x[], double y[], double *f, double norm )
+{
+//    x[0] -= y[0];
+//    x[1] -= y[1];
+    x[0] += y[0];
+    x[1] += y[1];
+    *f = Function(x[0], x[1]);
+    return 1;
+}
+
+
 double Step_subdivide( double x[], double y[], double *f, double norm )
 {
     double step = 1, fnew, epsilon = eps;
 }
 
 
-double Step_Newton( double x[], double y[], double *f, double norm )
-{
-    x[0] -= y[0];
-    x[1] -= y[1];
-    *f = Function(x[0], x[1]);
-    return 1;
-}
-
-
 //double lipsh = 9.2;
 
 //double Step_fixed( double x[], double y[], double *f, double norm )
 int Test( double (*Step)(double[], double[], double *, double), int i )
 {
     int    step = 0;
-    double x[2] = {-1, -1}, x_prev[2], delta[2], old_delta[2], y[2],
+    double x[2] = {/*0.2, 0.3*//*-0.5, 0.6*/1, 1}, x_prev[2], delta[2], old_delta[2], y[2],
             f = Function(x[0], x[1]), norm, mod, prev_mod;
     static double solution[2];
+    double antiGradient[2], curStep, fnew;
     
+    // Newton method with subdivision
+    if (i == 3)
+    {
+//        Gradient(x, y);
+//        mul(-1.0, y, antiGradient);
+//        curStep = chooseStepBySubdivision(x, &f, antiGradient, 0.5);
+//        printf("Initial step = %6.4lf\n\n", curStep);
+        
+        printf("x = (% 7.4lf, % 7.4lf) y = (% 7.4lf, % 7.4lf), f=%7.4lf \n",
+               x[0], x[1], y[0], y[1], f);
+        
+        do
+        {
+//            Gradient2(x, y, 1);
+                        
+            double H[2][2], invH[2][2], grad[2], det;
+            
+            Hessian(x, H);
+            det = H[0][0] * H[1][1] - H[0][1] * H[1][0];
+            if (det < 1e-5)
+            {
+                y[0] = y[1] = 0;
+                printf("det < 1e-5\n");
+            }
+            
+            invH[0][0] =  H[1][1] / det;
+            invH[0][1] = -H[0][1] / det;
+            invH[1][0] = -H[1][0] / det;
+            invH[1][1] =  H[0][0] / det;
+
+            Gradient(x, grad);
+            // Anti gradient
+            mul(-1.0, grad, grad);
+            
+            // Calc final direction
+            y[0] = (invH[0][0] * grad[0] + invH[0][1] * grad[1]);
+            y[1] = (invH[1][0] * grad[0] + invH[1][1] * grad[1]);
+            
+            curStep = chooseStepBySubdivision(x, &f, y, 0.5);
+            
+            mul(curStep, y, y);
+            Step(x, y, &f, norm);
+            
+            printf("x = (% 7.4lf, % 7.4lf) y = (% 7.4lf, % 7.4lf), f=%7.4lf ",
+                   x[0], x[1], y[0], y[1], f);
+            printf("step = %6.4lf\n", curStep);
+            
+            step++;
+            
+            fnew = Function(x[0] + y[0], x[1] + y[1]);
+        } while (fabs(fnew - f) > eps && step < 15); 
+        //while abs(f(vectorSum(x, vectorMultiply(direction, step))) - f(x)) > precision:
+        
+    }
+    else {
     do
     {
+        
         Gradient(x, y);
-        if (i == 3)
-            Gradient2(x, y, 1);
-//        if (i == 4)
-//            Gradient2(x, y, (step % period) == 0);
+        
+        //        if (i == 4)
+        //            Gradient2(x, y, (step % period) == 0);
         norm = sqrt(y[0] * y[0] + y[1] * y[1]);
         printf("x = (% 7.4lf, % 7.4lf) y = (% 7.4lf, % 7.4lf) norm = % 7.4lf,"
                " f=%7.4lf ", x[0], x[1], y[0], y[1], norm, f);
             printf(", delta = (%8.5lf, %8.5lf)", delta[0], delta[1]);
             if (step != 0)
                 printf(", delta * old_delta = %.5lf", delta[0] * old_delta[0] +
-                        delta[1] * old_delta[1]);
+                       delta[1] * old_delta[1]);
         }
         printf("\n");
-        
-        step++;
+
         old_delta[0] = delta[0];
         old_delta[1] = delta[1];
         prev_mod = mod;
+        
+        step++;
+        
     } while (norm > eps);
+    }
     
     if (i == 0)
         memcpy(solution, x, 2 * sizeof(double));
     printf("\n2-dimensional minimization example\n");
 //    printf("\nFixed step:\n");
 //    Test(Step_fixed, 0);
+    
+//    printf("TEST: %f\n\n", pow(2.0, 10));
+    
     printf("\nStep subdivision:\n");
     Test(Step_subdivide, 1);
     printf("\nOptimal step:\n");
     
     printf("\nSecond order gradient:\n");
     Test(Step_Newton, 3);
+    
 //    period = 2;
 //    for (int i = 0; i < 4; i++)
 //    {