如何用GSL求解变量数多于方程数的矩形非线性方程组？

问题背景 #

我用C++编写代码求解非线性方程组，已通过GSL实现了方程数m等于变量数n（m=n）的情况，但现在经常遇到方程数少于变量数（m<n）的矩形方程组，有以下疑问：

1. 能否直接用GSL现有功能求解这类矩形方程组？
2. 如果不行，必须采用其他方法或更换库吗？
3. 为什么GSL官方文档里没有关于求解矩形方程组的章节？

尝试的代码 #

#include <cmath>
#include <gsl/gsl_vector.h>
#include <iostream>
#include <gsl/gsl_multiroots.h>

void print_vector(gsl_vector *gsl) {
  printf("[");
  for (int r = 0; r < gsl->size; r++) {
    if (r != gsl->size - 1) {
      printf("%+.5f; ", gsl->data[r]);
    } else {
      printf("%+.5f", gsl->data[r]);
    }
  }
  printf("]\n");
}

int evaluate(const gsl_vector *x, void *params, gsl_vector *f) {
  const double x0 = gsl_vector_get(x, 0);
  const double x1 = gsl_vector_get(x, 1);
  const double x2 = gsl_vector_get(x, 2);

  const double y0 = 0.1 * pow(x0, 2) + 0.1 * pow(x1, 2) + 2 - x2;

  // 单独设置这一行无法运行
  gsl_vector_set(f, 0, y0);

  // 加上下面两行也无法运行
  // gsl_vector_set(f, 1, 0);
  // gsl_vector_set(f, 2, 0);

  // 加上下面两行同样无法运行
  // gsl_vector_set(f, 1, y0);
  // gsl_vector_set(f, 2, y0);

  return GSL_SUCCESS;
}

int main() {
  const size_t n = 3;

  const gsl_multiroot_fsolver_type *T;
  T = gsl_multiroot_fsolver_hybrids; // 报错："iteration is not making progress towards solution"
//  T = gsl_multiroot_fsolver_hybrid; // 同样报错迭代无进展
//  T = gsl_multiroot_fsolver_broyden; // 报错："gsl: lu.c:266: ERROR: matrix is singular"
//  T = gsl_multiroot_fsolver_dnewton; // 报错："gsl: lu.c:147: ERROR: matrix is singular"
  gsl_multiroot_fsolver *s = gsl_multiroot_fsolver_alloc(T, n);

  gsl_multiroot_function f_multiroot_function = {evaluate, n, nullptr};

  gsl_vector *x_gsl = gsl_vector_alloc(n);
  gsl_vector_set(x_gsl, 0, 3.4);
  gsl_vector_set(x_gsl, 1, 2.7);
  gsl_vector_set(x_gsl, 2, 3.8);

  std::cout << "x_gsl: ";
  print_vector(x_gsl);

  gsl_multiroot_fsolver_set(s, &f_multiroot_function, x_gsl);

  int status;
  size_t iter = 0;
  do {
    iter++;
    status = gsl_multiroot_fsolver_iterate(s);

    std::cout << "s->x=" << std::endl;
    print_vector(s->x);

    if (status) {
      std::cout << "STOPPING: " << gsl_strerror(status) << "." << std::endl;
      break;
    }

    status = gsl_multiroot_test_residual(s->f, 1e-5);
  } while (status == GSL_CONTINUE && iter < 20);

  return 0;
}

运行结果 #

结果1：迭代无进展 #

x_gsl: [+3.40000; +2.70000; +3.80000]
s->x=
[+3.40000; +2.70000; +3.80000]
s->x=
[+3.40000; +2.70000; +3.80000]
s->x=
[+3.40000; +2.70000; +3.80000]
s->x=
[+3.40000; +2.70000; +3.80000]
s->x=
[+3.40000; +2.70000; +3.80000]
s->x=
[+3.40000; +2.70000; +3.80000]
s->x=
[+3.40000; +2.70000; +3.80000]
s->x=
[+3.40000; +2.70000; +3.80000]
s->x=
[+3.40000; +2.70000; +3.80000]
STOPPING: iteration is not making progress towards solution.

Process finished with exit code 0

结果2：矩阵奇异报错 #

x_gsl: [+3.40000; +2.70000; +3.80000]
gsl: lu.c:147: ERROR: matrix is singular
Default GSL error handler invoked.

Process finished with exit code 134 (interrupted by signal 6: SIGABRT)

解答 #

1. GSL能否直接求解m<n的矩形非线性方程组？ #

不行。GSL的multiroots模块专为求解m=n的方形方程组设计，所有内置求解器（如你尝试的hybrids、broyden等）都依赖雅可比矩阵（或近似雅可比矩阵）可逆的前提——而m<n时，雅可比矩阵行不满秩，必然奇异，这就是你遇到"matrix is singular"错误的核心原因。

你尝试的补0、重复方程的做法没有意义：补0是引入无效约束，重复方程是冗余约束，都无法让求解器找到合理的搜索方向，只会导致迭代停滞或矩阵奇异。

2. 替代方案 #

m<n时，方程组的解不是唯一的，是一个(n-m)维的解空间。要得到具体解，有两种核心思路：

• 最小二乘路径：使用GSL的gsl_multifit_nlin非线性最小二乘模块，它专门处理m≠n的情况，本质是最小化残差的L2范数，找到最接近满足所有方程的解。
• 补充约束：根据实际业务逻辑，添加(n-m)个合理约束（比如变量取值范围、变量间的关系等），将矩形方程组转为方形方程组后，再用GSL的multiroots模块求解。
• 更换库：如果需要更灵活的求解能力，可以考虑Eigen、NLopt或Ipopt这类支持非线性最小二乘、约束优化的库。

3. 为什么GSL文档没有矩形方程组的章节？ #

GSL的模块划分逻辑清晰：

• multiroots负责求解有唯一精确根的方形方程组；
• multifit_nlin负责求解无唯一解、需在优化准则下找近似解的情况（含m<n的矩形方程组）。

矩形非线性方程组不存在"精确求解"的概念，只能在某种优化规则下寻找最优解，这属于最小二乘或约束优化范畴，而非根求解。GSL将这部分功能放在了非线性最小二乘模块，而非multiroots模块，因此你在multiroots文档中找不到相关内容。

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