我用最速下降法写了一个简单的梯度下降算法。
对于陡降，i意味着步长被优化为最小化f（x-lambda * grad（f））的步长，其中lambda是步长，结果是每个方向与前一个方向正交。
问题是这个程序看起来很不稳定：它只适用于二次函数，但即使是四次函数，如x^4 + y^4 + z^4，无论使用何种精度都是不稳定的

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

double l2_norm(double *x1, double *x2, int m);
void get_grad(double (*f)(double *x, int m), double *x, double *grad, int m);
void copy(double *a, double *b, int m);
void add(double *a, double *b, double lambda, int m);
double dFdLambda(double (*f)(double *x, int m), double *x, double *u, int m);
void descent(double (*f)(double *x, int m), double *x, double *u, int m);

void print_vec(double *x, int m){
    for(int i=0; i<m; i++){
        printf("%f\n", x[i]);
    }
    return;
}

double f_temp(double *x, int m){

    double f=0.;

    for(int i=0; i<m; i++){
        f += (x[i]-1)*(x[i]-1);
    }
    return f;
}

int main(){

    double x[] = {2., 2., 2.};
    double *grad;
    double *x_old;
    int m = 3;
    double e = 1.E-6;

    grad = malloc(m*sizeof(double));
    x_old = malloc(m*sizeof(double));
    
    while(l2_norm(x, x_old, m) > e){
        copy(x, x_old, m);
        get_grad(f_temp, x, grad, m);
        descent(f_temp, x, grad, m);
        print_vec(x, m);
    }

    printf("\n");
    print_vec(x, m);
    

    return 0;
}

double l2_norm(double *x1, double *x2, int m){

    double norm = 0;

    for(int i=0;i<m;i++){
        norm += pow(x1[i]-x2[i], 2);
    }

    norm = sqrt(norm);

    return norm;
}

void get_grad(double (*f)(double *x, int m), double *x, double *grad, int m){
    /*
    numerical gradient with simmetric method

    the actual gradient computed is -grad in order to be used in 
    gradient descent
    */
    double e = 1.E-6; // numerical precision
    double *x_forward;
    double *x_backward;
    x_forward = malloc(m*sizeof(double));
    x_backward = malloc(m*sizeof(double));
    
    for(int i=0;i<m; i++){
        x_forward[i] = x[i] + e;
        x_backward[i] = x[i] - e;    
        grad[i] = -(f(x_forward, m) - f(x_backward, m))/(2*e);
        x_forward[i] -= e;
        x_backward[i] += e;
    }
    free(x_forward);
    free(x_backward);
    return;
}

double dFdLambda(double (*f)(double *x, int m), double *x, double *u, int m){
    /*
    compute the total derivative dF/dLambda in the point x along
    the direction u
    */
    double e = 1.E-5;
    double *x_forward;
    double *x_backward;
    double der;

    x_forward = malloc(m*sizeof(double));
    x_backward = malloc(m*sizeof(double));
    copy(x, x_forward, m);
    copy(x, x_backward, m);
    add(x_forward, u, e, m);
    add(x_backward, u, -e, m);

    der = (f(x_forward, m) - f(x_backward, m))/(2*e);

    free(x_forward);
    free(x_backward);

    return der;
}

void copy(double *a, double *b, int m){
    //copies a into b
    for(int i=0; i<m; i++){
        b[i] = a[i];
    }
    return;
}

void add(double *a, double *b, double lambda, int m){
    /*
    adds lambda*b to a
    */
    for(int i=0; i<m; i++){
        a[i] += lambda*b[i];
    }
    return;
}

void descent(double (*f)(double *x, int m), double *x, double *u, int m){
    /*
    actual gradient descent starting from x going in direction u
    */
    double e = 1.E-5;
    double derA, derB, derC;
    double *x_start;
    double lambda, lambda_min, lambda_max;

    lambda = e;

    x_start = malloc(m*sizeof(double));
    copy(x, x_start, m);
    derA = dFdLambda(f, x, u, m);
    derC = derA;

    /*
    this while loop finds the interval in which the total derivatives df/dl
    changes sign i.e. the interval in which the solution of df/dl = 0
    will be searched with bisection method
    */
    while((derA*derC) >= 0.){
        copy(x_start, x, m);
        add(x, u, lambda, m);
        derC = dFdLambda(f, x, u, m);
        lambda *= 2.;
    }

    /*
    because lambda >= 0 the leftmost point of the interval is 0, the rightmost 
    the point found before after which the total derivatives changes sign
    */
    lambda_min = 0;
    lambda_max = lambda;
    lambda = 0.5*(lambda_min + lambda_max);

    while((fabs(lambda_max - lambda_min)) > e){
        copy(x_start, x, m);
        add(x, u, lambda_min, m);
        derA = dFdLambda(f, x, u, m);

        copy(x_start, x, m);
        add(x, u, lambda, m);
        derB = dFdLambda(f, x, u, m);

        if((derA*derB) > 0.){
            lambda_min = lambda;
        } 
        else if((derA*derB) < 0.){
            lambda_max = lambda;
        }
        else{
            lambda_min = lambda_max = lambda;
        }

        lambda = 0.5*(lambda_min + lambda_max);
    }

    copy(x_start, x, m);
    // the modified vector x will cointain the coordinates of the minumum
    add(x, u, lambda, m); 
    free(x_start);

    return;
}

使用f += (x[i]-1)*(x[i]-1)，输出正确

1.000005
1.000005
1.000005
-1.499995
-1.499995
1.000005
1.000012
1.000012
1.000000
1.000000
1.000000
1.000000
1.000000
1.000000
1.000000

1.000000
1.000000
1.000000

但是使用(x[i]-1)*(x[i]-1)*(x[i]-1)*(x[i]-1)或者甚至使用幂为2的pow（）函数，结果类似于

1.000010
1.000010
1.000010
-1.499990
-1.499990
1.000010
1.000292
1.000292
1.000010
1.000292
1.000000
1.000010
1.000292
-1032674302637538115769421495899732814941072310657452010951926707678661836800.000000
-171980112174886070524247400874370796081962558827236065865505533642887865040896.000000
inf
inf
nan

inf
inf
nan

降低函数（double e变量）内部的数值精度似乎只适用于简单的二次函数（不使用pow（）函数），否则我会得到null或inf数字