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

How to implement a left matrix division on C++ using gsl

本文关键字:除法 实现 何使用 gsl c++      更新时间:2023-10-16

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

Am-by-n矩阵,其中m不等于n, B是具有m分量的列向量。

并且我希望结果X = AB是欠定或过定方程组AX = B的最小二乘解。换句话说,X最小化了向量AX - B的长度norm(A*X - B)。这意味着我希望它具有与MATLAB中的AB相同的结果。

我想在GSL-GNU (GNU Science Library)中实现这个功能,我不太了解数学,最小二乘拟合或矩阵运算,有人能告诉我如何在GSL中做到这一点吗?或者如果在GSL中实现它们太复杂,有人能建议我一个好的开源C/c++库,提供上述矩阵操作吗?


好吧,我又花了5个小时才自己弄明白。但还是要感谢你对我问题的建议。

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

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

和向量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;
 }

除了您给出的例子,快速搜索显示了其他GSL示例,一个使用QR分解,另一个使用LU分解。

存在其他能够求解线性系统的数值库(每个线性代数库的基本功能)。例如,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是一个超级函数,具有多个执行路径。如果系数矩阵A具有某种特殊的结构,则利用它可以获得更快或更准确的结果(可以选择代换算法、LU和QR分解等)

MATLAB也有PINV,它返回最小范数最小二乘解,除了许多其他迭代方法用于求解线性方程组。

我不确定我是否理解你的问题,但如果你已经找到了使用MATLAB的解决方案,你可能要考虑使用MATLAB编码器,它会自动将你的MATLAB代码翻译成c++。