C语言实现梯度下降算法的数值不稳定问题咨询

我实现了带步长优化的最速下降法（通过优化步长λ使得f(x - λ·∇f)取最小值，保证迭代方向正交），但程序仅在二次函数上能正常运行，对x⁴+y⁴+z⁴这类四次函数，无论调整精度都会出现数值发散，输出无穷大或NaN。

实现代码 #

#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)后，输出出现数值爆炸：

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）仅对简单二次函数有效，其他情况仍会出现NaN、无穷大等异常。求解决该数值不稳定问题的方案。

内容的提问来源于stack exchange，提问作者Niccolò Tiezzi