C中梯度下降的数值不稳定性

y53ybaqx  于 2023-06-21  发布在  其他
关注(0)|答案(2)|浏览(99)

我用最速下降法写了一个简单的梯度下降算法。
对于陡降,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数字

c7rzv4ha

c7rzv4ha1#

在使用值不确定的已分配对象时,您会遇到一些问题。第一个是main

x_old = malloc(m*sizeof(double));
    
    while(l2_norm(x, x_old, m) > e){

请注意,您使用的是x_old指向的数据,而没有为其分配任何值。由于您总是希望至少执行一次迭代,因此将while循环重构为do ... while循环是有意义的。
但是你在get_grad()中还有一些可能更有影响力的东西:

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;
    }

请注意,当x_forwardx_backward的值未定义时,在分配后立即开始使用这些值所指向的数据。你 * 最终 * 为所有分配的元素赋值,但你 * 最初 * 只为每个元素赋值0。假设您希望将x中的数据复制到开始循环之前这些数据所指向的空间中。如目前所写的,在没有这样的复制的情况下,所得到的行为是未定义的。绝对没有理由期望存储在*grad中的结果数据将是x处梯度的估计值。
即使您可以依靠malloc()将分配的空间初始化为全零(它不这样做,但calloc()这样做),程序也会出错。特别是,get_grad()函数通常会计算假设梯度在不同点的每个分量。
我能相当精确地重现你程序的错误输出。计算采取了类似的路径,在不太多的迭代之后发散到NaN。
在解决了这两个问题之后,程序会生成以下输出:

1.000005
1.000005
1.000005
1.000000
1.000000
1.000000
1.000000
1.000000
1.000000

1.000000
1.000000
1.000000

更新

该程序在函数dFdLambda()中也有一个实现错误。它计算 f 在所选方向上的导数为

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

,但是那个公式中的除数是不正确的,在某些情况下严重错误。它应该是x_forwardx_backward之间的距离(L2范数)。当我做了这个修正后,得到的程序在四次情况下也收敛了。
我想你可能仍然可以找到输入函数的程序行为不当。特别是,对于在局部最小值附近具有更陡梯度的函数,以及那些振荡非常快的函数,您可能会遇到麻烦。对于在域的各个维度上具有明显各向异性的函数,以及具有非常大的值的局部极小值的函数,您也可能会遇到麻烦。
解决@chux描述的问题将有助于解决其中的一些问题。一个更适应的方法来估计梯度和导数可能会帮助一些其他的。
下降计算本身可以防止发散。例如,可以验证 f 的值实际上确实减小,和/或可以以各种方式抑制或限制每个下降步长的 x 的变化。
此外,我理解 * 数学 * 为什么descent()估计函数的导数并使用它们沿着所选方向搜索最小值,通过 * 计算 *,您可能会更好地直接搜索最小值,使用函数的值而不是其导数的估计值。这肯定会减少你所观察到的分歧的可能性,而且它的收敛特性一般来说可能一样好或更好。
然而,最终,数值程序受到它们使用的数值表示的限制,这些数值表示不可避免地具有有限的范围和精度(尽管如果你愿意,你可以或多或少地任意选择特定的限制)。这是数值编程的一个基本约束。

  • 虽然我不确定为什么它会有区别,因为程序似乎只依赖于导数的符号,而不是它的大小。
wwwo4jvm

wwwo4jvm2#

除了@John Bollinger很好的回答:

考虑一个替代的“亲密度”测试

无论使用何种精度都不稳定
代码目前比较线性意义上的值,但浮点值是对数分布的-这就是为什么它们被称为 floating point而不是 fixed point。他们应该得到一个使用根据参数大小而变化的测试的测试。否则,所有的大数对永远不会“大约相等”,除非完全相等,而所有的小数对总是相等的。

double e = 1.E-5;
...
// while((fabs(lambda_max - lambda_min)) > e){
while(!about_equal(lambda_max, lambda_min, e)) {

下面的候选about_equal()是 * 说明性的 *,没有优化。如何进行最佳比较取决于OP目标的细节。

bool about_equal(double a, double b, double relative_error) {
  double diff = fabs(a - b);
  double magnitude = fmax(fabs(a), fabs(b));
  return diff <= magnitude * relative_error;
}

相关问题