Conjugate Gradient Method + Strong Wolfe Condition

CGSW.c

#include <stdio.h>
#include <math.h>
#define N 21                    // number of variables +1 (shouble be odd)

double delta_1 = 1E-4;          // strong wolfe condition Rule1 delta_1
double delta_2 = 0.1;           // strong wolfe condition Rule2 delta_2
double epsilon = 1E-6;

// define f function
double f(double x[N]) {
    double y = 0;
    int i = 1;
    for (i = 1; i <= 10; i++) {
        y += (1 - x[2 * i - 1]) * (1 - x[2 * i - 1]) + 10 * (x[2 * i] - x[2 * i - 1] * x[2 * i - 1])
             * (x[2 * i] - x[2 * i - 1] * x[2 * i - 1]);
    }
    return y;
}

// f function's gradient
void Gradient(double x[N], double g_f[N]) {
    int i = 1;
    for (i = 1; i <= 10; i++) {
        g_f[2 * i - 1] = 2 * (x[2 * i - 1] - 1) + 40 * (x[2 * i - 1] * x[2 * i - 1] - x[2 * i])
                         * x[2 * i - 1];
        g_f[2 * i] = 20 * (x[2 * i] - x[2 * i - 1] * x[2 * i - 1]);
    }
}

// gradient's norm
double Gradient_f_norm(double g_f[N]) {
    double gradient_f_norm = 0;
    int i;
    for (i = 1; i < N; i++) {
        gradient_f_norm += g_f[i] * g_f[i];
    }
    return sqrt(gradient_f_norm);

}

// beta_k(FR)
// g_f1 stands for the one before g_f => (x1 stands for the one before x)
double Beta_k(double g_f_norm, double g_f1_norm) {
    return g_f_norm * g_f_norm / (g_f1_norm * g_f1_norm);
}

// direction:d_k
void D_k(double *d_k, double *d_k1, int k, double *x_k, double *x_k1, double g_f[N], double beta_k) {
    int i;
    if (k == 0)
        for (i = 1; i < N; i++) {
            d_k[i] = -g_f[i];
            d_k1[i] = d_k[i];
        }
    else
        for (i = 1; i < N; i++) {
            d_k[i] = -g_f[i] + beta_k * d_k1[i];
            d_k1[i] = d_k[i];
        }
}

// update x_k:x_(k+1)=x_k+alpha_k * d_k
void X_k_update(double *x_k, double *x_k1, double alpha_k, double *d_k) {
    int i = 1;
    for (i = 1; i < N; i++) {
        x_k1[i] = x_k[i];
        x_k[i] = x_k[i] + alpha_k * d_k[i];
    }
}

// dot product of g_f(gradient) and d_k(direction)
double Dotproduct_gradient_d_k(double x_k[N], double d_k[N]) {
    double dotproduct = 0.0, g_f[N];
    Gradient(x_k, g_f);
    int i;
    for (i = 1; i < N; i++) {
        dotproduct += g_f[i] * d_k[i];
    }
    return dotproduct;
}

// strong wolfe condition rule1
double Rule_1(double *x_k, double *d_k, double alpha_k) {
    double x_k_temp[N];
    int i;
    for (i = 1; i < N; i++) {
        x_k_temp[i] = x_k[i] + alpha_k * d_k[i];
    }
    if (f(x_k_temp) <= f(x_k) + delta_1 * alpha_k * Dotproduct_gradient_d_k(x_k, d_k))
        return 1;
    else return 0;
}

// strong wolfe condition rule2
double Rule_2(double *x_k, double *d_k, double alpha_k) {
    double x_k_temp[N];
    int i;
    for (i = 1; i < N; i++) {
        x_k_temp[i] = x_k[i] + alpha_k * d_k[i];
    }
    if (fabs(Dotproduct_gradient_d_k(x_k_temp, d_k)) <= delta_2 * fabs(Dotproduct_gradient_d_k(x_k, d_k)))
        return 1;
    else return 0;
}

// stepsize:search for alpha_k(fit strong wolfe condition)
double Strong_wolfe(double *x_k, double *d_k, double g_f[N]) {
    double alpha = 1, alpha_min = 0, alpha_max = 1, alpha_limit = 0.1;
    double ak = alpha_min;
    double bk = alpha_max;
    double alphap;
    double fptemp;
    double f1temp = f(x_k);
    double ftemp;
    double f1ptemp = Dotproduct_gradient_d_k(x_k, d_k);
    double x_k_temp[N];
    double bracketSize;
    int i;
    while (1) {
        for (i = 1; i < N; i++) {
            x_k_temp[i] = x_k[i] + alpha * d_k[i];
        }
        ftemp = f(x_k_temp);
        if (Rule_1(x_k, d_k, alpha)) {
            for (i = 1; i < N; i++) {
                x_k_temp[i] = x_k[i] + alpha * d_k[i];
            }
            fptemp = Dotproduct_gradient_d_k(x_k_temp, d_k);
            if (Rule_2(x_k, d_k, alpha)) {
                //                printf("Yes! Found one alpha_k!\n");
                return alpha;
            } else {
                alphap = alpha + (alpha - ak) * fptemp / (f1ptemp - fptemp);
                ak = alpha;
                alpha = alphap;
                f1temp = ftemp;
                f1ptemp = fptemp;
            }
        } else {
            bracketSize = fabs(bk - ak);
            alphap = ak + (alpha - ak) / (2 * (1 + (f1temp - ftemp) / ((alpha - ak) * f1ptemp)));
            bk = alpha;
            if (fabs(alphap - ak) < alpha_limit * bracketSize) {
                alpha = ak + alpha_limit * bracketSize;
            } else {
                alpha = alphap;
            }
        }
    }

}

// Array x[N] <= y[N]
void Trans_array(double x[N], double y[N]) {
    int i;
    for (i = 1; i < N; i++) {
        x[i] = y[i];
    }
}

// print x_k
void Print_xk(double *x_k) {
    int i;
    printf("x_k:\n");
    for (i = 1; i < N; i++) {
        printf("x_k[%d]:%6.4lf", i, x_k[i]);
        if (i % 4 == 0)
            printf("\n");
        else
            printf("\t");
    }
}

// print more accurate x_k
void Print_xk_accurate(double *x_k) {
    int i;
    printf("x_k:\n");
    for (i = 1; i < N; i++) {
        printf("x_k[%d]:%6.9lf", i, x_k[i]);
        if (i % 2 == 0)
            printf("\n");
        else
            printf("\t");
    }
}

// print key variables & parameters
void Print_state(double x_k[N], double g_f_norm, double beta_k, double alpha_k, int k) {
    printf("===============================================================\n");
    printf("k=%d\n", k);
    printf("-----\n");
    Print_xk(x_k);
    printf("-------------------------------------------------------------\n");
    printf("alpha_k:%lf      beta_k:%lf      g_f_norm:%lf\n", alpha_k, beta_k, g_f_norm);
    printf("-------------------------------------------------------------\n");
    printf("f(x_k) = %0.15lf\n", f(x_k));
}

int main() {
    int i, k = 0;
    double x_k[N];
    // double x_k[N] = {1, -1.2, 1, -1.2, 1, -1.2, 1, -1.2, 1, -1.2, 1, -1.2, 1, -1.2, 1, -1.2, 1, -1.2, 1, -1.2, 1};
    double x_k1[N], d_k[N], d_k1[N], g_f[N], g_f1[N], g_f_norm, g_f1_norm, alpha_k = 1.0, beta_k = 1;
    printf("Please type in the initial value of these variables......\n");
    for (i = 1; i < N; i++) {
        scanf("%lf", &x_k[i]);
    }
    printf("initial value of x:\n");
    Print_xk(x_k);
    printf("-------------------------------------------------------------\n");
    printf("f(x_k)=%0.15lf\n", f(x_k));
    Gradient(x_k, g_f);
    D_k(d_k, d_k1, k, x_k, x_k1, g_f, beta_k);                          // get first d_k
    k = k + 1;
    Trans_array(x_k1, x_k);
    Trans_array(g_f1, g_f);
    g_f_norm = Gradient_f_norm(g_f);
    while (g_f_norm > epsilon) {
        alpha_k = Strong_wolfe(x_k, d_k, g_f);                          // get alpha_k with strong wolfe
        X_k_update(x_k, x_k1, alpha_k, d_k);                            // x_k_update
        Gradient(x_k, g_f);
        g_f1_norm = g_f_norm;
        g_f_norm = Gradient_f_norm(g_f);
        beta_k = Beta_k(g_f_norm, g_f1_norm);                           // get beta_k
        D_k(d_k, d_k1, k, x_k, x_k1, g_f, beta_k);                      // get next d_k
        Print_state(x_k, g_f_norm, beta_k, alpha_k, k);
        k += 1;
    }
    printf("*************************************************************\n");
    printf("Finally, x_star:\n");
    printf("------------------\n");
    Print_xk_accurate(x_k);
    printf("-------------------------------------------\n");
    printf("g_f_norm = %8.7lf < epsilon = %8.7lf\n", g_f_norm, epsilon);
    printf("-------------------------------------------\n");
    printf("f(x_k) = %0.15lf, k = %d\n", f(x_k), k - 1);
    printf("*************************************************************\n");
    return 0;
}