Julia中多变量牛顿法调用时SingularException(2)问题排查

算力阿峰

2026-5-12

多变量牛顿法遇到SingularException及异常结果的原因与解决办法

问题复现

你在Julia中实现了多变量牛顿法：

function newton(f::Function, J::Function, x)
    h = Inf64
    tolerance = 10^(-10)
    while (norm(h) > tolerance)
        h = J(x)\f(x)
        x = x - h
    end
    return x
end

当求解第一个方程组时抛出LoadError: SingularException(2)：

f(θ) = [cos(θ[1]) + cos(θ[2] - 1.3), sin(θ[1]) + sin(θ[2]) - 1.3]
J(θ) = [-sin(θ[1]) -sin(θ[2]); cos(θ[1]) cos(θ[2])]
θ = [pi/3, pi/7]
newton(f, J, θ)

而第二个方程组能正常返回正确解：

f(x) = [(93-x[1])^2 + (63-x[2])^2 - 55.1^2, (6-x[1])^2 + (16-x[2])^2 - 46.2^2]
J(x) = [-2*(93-x[1]) -2*(63-x[2]); -2*(6-x[1]) -2*(16-x[2])]
x = [35, 50]
newton(f, J, x)

另外，先求解第二个方程组再求解第一个时，前者不再报错但结果完全不准确。

原因分析

1. 第一个方程组的SingularException：雅可比矩阵定义错误

你的第一个方程组的雅可比矩阵计算有误：
第一个方程是 cos(θ[1]) + cos(θ[2] - 1.3)，对θ[2]的偏导数应该是 -sin(θ[2] - 1.3)，但你写成了 -sin(θ[2])。错误的雅可比矩阵会导致迭代过程中，某一步的矩阵行列式为0（奇异矩阵），从而无法进行线性方程组求解（J(x)\f(x)本质是求逆运算），抛出SingularException。
即使初始点的雅可比矩阵行列式不为0，错误的偏导数也会让牛顿法的迭代方向完全偏离正确路径，最终大概率会遇到奇异矩阵的情况。

2. 先后运行后的异常结果：函数名覆盖导致调用错误

Julia中，后定义的同名函数会覆盖之前的定义。如果你在同一个会话中先运行第二个方程组的代码（定义了f(x)和J(x)），再运行第一个方程组时没有重新定义f和J，直接调用newton(f, J, θ)，此时实际调用的是第二个方程组的f和J函数，而非第一个的。这就导致你在求解第一个问题时，实际计算的是第二个方程组在θ=[pi/3, pi/7]处的迭代结果，自然返回的解完全不准确，且因为第二个的雅可比矩阵在该点可逆，所以不会抛出奇异错误。

解决办法

1. 修正第一个方程组的雅可比矩阵

将J(θ)的第一行第二列修正为-sin(θ[2] - 1.3)，正确的雅可比矩阵代码如下：

J(θ) = [-sin(θ[1])  -sin(θ[2] - 1.3); 
         cos(θ[1])   cos(θ[2])]

修正后，重新定义f和J并调用牛顿法，就能正常迭代收敛到正确解。

2. 避免函数名冲突（可选）

为了避免不同方程组的函数名覆盖问题，可以给不同的方程组定义不同的函数名，比如第一个用f1、J1，第二个用f2、J2，这样即使在同一个会话中先后运行，也不会出现调用错误：

# 第一个方程组
f1(θ) = [cos(θ[1]) + cos(θ[2] - 1.3), sin(θ[1]) + sin(θ[2]) - 1.3]
J1(θ) = [-sin(θ[1])  -sin(θ[2] - 1.3); 
         cos(θ[1])   cos(θ[2])]
θ = [pi/3, pi/7]
newton(f1, J1, θ)

# 第二个方程组
f2(x) = [(93-x[1])^2 + (63-x[2])^2 - 55.1^2, (6-x[1])^2 + (16-x[2])^2 - 46.2^2]
J2(x) = [-2*(93-x[1]) -2*(63-x[2]); -2*(6-x[1]) -2*(16-x[2])]
x = [35, 50]
newton(f2, J2, x)

3. 添加迭代保护（可选）

为了让牛顿法更健壮，可以在迭代过程中检查雅可比矩阵的奇异性，比如添加行列式判断或使用伪逆替代直接求逆：

function robust_newton(f::Function, J::Function, x)
    h = Inf64
    tolerance = 10^(-10)
    max_iter = 1000
    iter = 0
    while (norm(h) > tolerance) && (iter < max_iter)
        jac = J(x)
        # 检查行列式是否接近0，避免奇异矩阵
        if abs(det(jac)) < 1e-8
            println("Warning: Near-singular Jacobian at iteration $iter, using pseudoinverse")
            h = pinv(jac) * f(x)
        else
            h = jac \ f(x)
        end
        x = x - h
        iter += 1
    end
    if iter == max_iter
        println("Warning: Reached maximum iterations")
    end
    return x
end

这个版本的牛顿法在遇到接近奇异的矩阵时，会使用伪逆(pinv)继续迭代，避免直接抛出异常。

验证结果

修正雅可比矩阵后，调用newton(f, J, θ)，会收敛到正确的解（代入解到f(θ)中，结果会接近0向量）。

内容的提问来源于stack exchange，提问作者K. Claesson