c++ - 如何使用gsl在C++上实现左矩阵除法

Question

我正在尝试将 MATLAB 程序移植到 C++。 我想在矩阵 A 和列 vector B 之间实现左矩阵除法。

A 是一个 m-by-n 矩阵，其中 m 不等于 n 并且 B 是一个包含m 个分量的列 vector 。

我希望结果 X = A\B 是欠定或超定方程组 AX = B 的最小二乘解。换句话说，X 最小化 norm(A*X - B)，即 vector AX - B 的长度。 这意味着我希望它与 MATLAB 中的 A\B 具有相同的结果。

我想在 GSL-GNU(GNU 科学图书馆)中实现这个功能，但我不太了解数学、最小二乘拟合或矩阵运算，有人能告诉我如何在 GSL 中实现吗？或者，如果在 GSL 中实现它们太复杂，有人可以向我推荐一个提供上述矩阵运算的良好开源 C/C++ 库吗？

---

好吧，我又花了 5 个小时终于自己弄明白了。但仍然感谢您对我的问题提出的建议。

假设我们有一个 5 * 2 的矩阵

A = [1 0           
     1 0
     0 1
     1 1
     1 1] 

和一个 vector b = [1.8388,2.5595,0.0462,2.1410,0.6750]

A\b 的解决方案是

 #include <stdio.h>
 #include <gsl/gsl_linalg.h>

 int
 main (void)
 {
   double a_data[] = {1.0, 0.0,1.0, 0.0, 0.0,1.0,1.0,1.0,1.0,1.0};

   double b_data[] = {1.8388,2.5595,0.0462,2.1410,0.6750};

   gsl_matrix_view m
     = gsl_matrix_view_array (a_data, 5, 2);

   gsl_vector_view b
     = gsl_vector_view_array (b_data, 5);

   gsl_vector *x = gsl_vector_alloc (2); // size equal to n
   gsl_vector *residual = gsl_vector_alloc (5); // size equal to m
   gsl_vector *tau = gsl_vector_alloc (2); //size equal to min(m,n)
   gsl_linalg_QR_decomp (&m.matrix, tau); // 
   gsl_linalg_QR_lssolve(&m.matrix, tau, &b.vector, x, residual);

   printf ("x = \n");
   gsl_vector_fprintf (stdout, x, "%g");
   gsl_vector_free (x);
   gsl_vector_free (tau);
   gsl_vector_free (residual);
   return 0;
 }

Answer 1

除了你给的那个，快速搜索一下发现了其他GSL examples , 一个使用 QR decomposition , 其他 LU decomposition .

还有其他numeric libraries能够求解线性系统(每个线性代数库中的基本功能)。对于一个，Armadillo提供了一个漂亮且可读的界面:

#include <iostream>
#include <armadillo>
using namespace std;
using namespace arma;

int main()
{
    mat A = randu<mat>(5,2);
    vec b = randu<vec>(5);

    vec x = solve(A, b);
    cout << x << endl;

    return 0;
}

另一个不错的是 Eigen图书馆:

#include <iostream>
#include <Eigen/Dense>
using namespace std;
using namespace Eigen;

int main()
{
    Matrix3f A;
    Vector3f b;
    A << 1,2,3,  4,5,6,  7,8,10;
    b << 3, 3, 4;

    Vector3f x = A.colPivHouseholderQr().solve(b);
    cout << "The solution is:\n" << x << endl;

    return 0;
}

---

现在，要记住的一件事是 MLDIVIDE是 super-charged功能，并有多个执行路径。如果系数矩阵A有一些特殊的结构，那么利用它来获得更快或更准确的结果(可以选择替代算法，LU和QR分解，..)

MATLAB 也有 PINV它返回最小范数最小二乘解，以及其他一些 iterative methods用于求解线性方程组。