Julia梯度下降脚本多次执行结果异常：从正确值变为NaN

问题背景 #

尝试用Julia实现梯度下降优化算法，首次执行可得到正确结果（x=[-1.0,1.5]，α=1.0），但重复执行脚本时结果逐渐变化，最终全部变为NaN。推测是变量残留了上一次执行的值，但无法定位具体变量。原代码如下：

using LinearAlgebra #pacote de álgebra linear

# Parâmetros importantes do método do gradiente 
h = 0.001 #passo da derivada numérica
precision_a = 1.0e-04 #precisão da função de otimização de alfa
precision_fun = 1.0e-08 #precisão da função gradiente
a0 = 0.0 #chute inicial de alfa

#parâmetros de entrada
f1_jac(x1::Float64,x2::Float64) = 1 + 4x1 +2x2
f2_jac(x1::Float64,x2::Float64) = -1 + 2x1 +2x2
jac = [f1_jac,f2_jac]
n = length(jac)
x0 = [0.0,0.0]

#variáveis auxiliares
grad_a = zeros(n) 
grad_xk = zeros(n)
grad_xk1 = zeros(n)

function f_a(a::Float64,x::Array)
    for i=1:n
        grad_a[i] = jac[i](x...)
    end
    xk1 = x - a.*grad_a
    for i=1:n
        grad_xk1[i] = jac[i](xk1...)
    end
    return transpose(grad_a)*grad_xk1
end

function df_a(a::Float64,x::Array)
    df_a = (f_a(a+h,x)-f_a(a-h,x))/2h
    return df_a
end

function teste(x::Array)
    xa = x
    k_alfa = 0
    a_k = 0
    a_dummy = a0
    erro_a = 1.0
    while erro_a>precision_a && k_alfa<=1000
        k_alfa += 1
        a_k = a_dummy - f_a(a_dummy,xa)/df_a(a_dummy,xa)
        erro_a = abs(a_k-a_dummy)
    end
    return a_k
end

function teste2()
    k_it = 0
    x_dummy = x0 
    erro_fun = 1.0
    x_k = zeros(n)
    while erro_fun>precision_fun && k_it<=1000
        alfa = teste(x_dummy)
        x_k = x_dummy - alfa.*grad_a
        erro_fun = max(abs.(grad_a)...)
        x_dummy .= x_k
        k_it+=1
        println("$k_it = iteração e $(trunc(alfa,digits = 2)) = alfa")
    end
    return trunc.(x_k,digits=2)
end    

问题定位 #

1. 全局变量状态污染

   • grad_a、grad_xk1是全局数组，f_a函数直接修改这些变量的内容，每次执行后变量不会自动重置，残留的旧值会干扰后续计算。
   • x0作为全局初始值数组，teste2中x_dummy = x0是引用赋值，后续x_dummy .= x_k会直接修改x0的内容，导致下一次执行teste2时初始值不再是[0.0,0.0]，这是重复执行结果变化的核心原因。

2. 步长优化迭代逻辑错误

   • teste函数的while循环中，计算得到新的步长a_k后，未更新a_dummy的值，导致循环始终使用初始的a0=0.0进行迭代，无法推进步长优化，最终返回错误的α值。

3. 数值微分的潜在问题

   • df_a中/2h的写法容易引发运算优先级混淆（应明确写为/(2*h)），且固定步长h=0.001可能导致数值微分精度不足，影响步长优化的稳定性。

代码优化建议与修正后的代码 #

优化要点 #

• 将辅助变量移至函数内部，避免全局状态污染；
• 对初始值数组进行复制，保护全局初始值不被修改；
• 修复步长优化的迭代逻辑，更新迭代变量；
• 明确数值微分的运算优先级，提升代码可读性；
• 添加精确的类型标注，增强代码稳定性与性能。

修正后的代码 #

using LinearAlgebra

# 梯度下降核心参数
h = 0.001                  # 数值微分步长
precision_a = 1.0e-04      # 步长优化精度
precision_fun = 1.0e-08    # 梯度收敛精度
a0 = 0.0                   # 步长初始猜测值

# 梯度分量函数
f1_jac(x1::Float64, x2::Float64) = 1 + 4x1 + 2x2
f2_jac(x1::Float64, x2::Float64) = -1 + 2x1 + 2x2
jac = [f1_jac, f2_jac]
n = length(jac)
x0 = [0.0, 0.0]            # 全局初始值，仅作为模板使用

function f_a(a::Float64, x::Vector{Float64})
    # 局部初始化梯度变量，避免全局状态污染
    grad_a = zeros(Float64, n)
    grad_xk1 = zeros(Float64, n)
    
    for i in 1:n
        grad_a[i] = jac[i](x...)
    end
    xk1 = x - a .* grad_a
    for i in 1:n
        grad_xk1[i] = jac[i](xk1...)
    end
    return transpose(grad_a) * grad_xk1
end

function df_a(a::Float64, x::Vector{Float64})
    # 明确运算优先级，避免歧义
    return (f_a(a + h, x) - f_a(a - h, x)) / (2 * h)
end

function teste(x::Vector{Float64})
    xa = copy(x)  # 复制输入，避免修改原数组
    k_alfa = 0
    a_dummy = a0
    erro_a = 1.0
    
    while erro_a > precision_a && k_alfa <= 1000
        k_alfa += 1
        fa_val = f_a(a_dummy, xa)
        dfa_val = df_a(a_dummy, xa)
        # 避免除以0的情况
        if isapprox(dfa_val, 0.0; atol=1e-10)
            break
        end
        a_k = a_dummy - fa_val / dfa_val
        erro_a = abs(a_k - a_dummy)
        a_dummy = a_k  # 更新迭代变量，推进优化
    end
    return a_dummy
end

function teste2()
    k_it = 0
    x_dummy = copy(x0)  # 复制初始值，保护全局x0不被修改
    erro_fun = 1.0
    x_k = zeros(Float64, n)
    
    while erro_fun > precision_fun && k_it <= 1000
        k_it += 1
        alfa = teste(x_dummy)
        
        # 重新计算当前点的梯度（避免依赖全局变量）
        grad_current = zeros(Float64, n)
        for i in 1:n
            grad_current[i] = jac[i](x_dummy...)
        end
        
        x_k = x_dummy - alfa .* grad_current
        erro_fun = maximum(abs.(grad_current))
        x_dummy .= x_k
        
        println("$k_it = iteração e $(trunc(alfa, digits=2)) = alfa")
    end
    return trunc.(x_k, digits=2)
end

验证说明 #

修正后的代码每次执行teste2()时，都会从初始值[0.0,0.0]开始计算，所有辅助变量都在函数内部初始化，不会残留上一次的状态。步长优化的迭代逻辑正常推进，能够稳定得到正确结果[-1.0, 1.5]。

内容的提问来源于stack exchange，提问作者elrd12