How to implement splitting left matrix in C ++ using gsl

I am trying to port a MATLAB program to C ++. And I want to implement left matrix division between matrix Aand column vector B.

A- this is a matrix m-by-nwith mnot equal to n, but Bis a column vector with components m.

And I want the result to X = A\Bbe a least squares solution to an under- or overdetermined system of equations AX = B. In other words, Xminimizes the norm(A*X - B)length of the vector AX - B. This means that I want it to have the same result as A\Bin MATLAB.

I want to implement this function in the GSL-GNU (GNU Science Library), and I don’t know too much about math, least squares or matrix operation, can someone tell me how to do this in GSL? Or if it’s too difficult to implement them in GSL, can anyone suggest me a good open source C / C ++ library that provides matrix operation?


Well, I finally found out that I spent another 5 hours on this. But thanks for the suggestions on my question.

Assuming we have a 5 * 2 matrix

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

and vector b = [1.8388,2.5595,0.0462,2.1410,0.6750]

Decision A \ bwill be

 #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;
 }
+5
source share
2 answers

, , 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, , .

+4

, , MATLAB, MATLAB Coder, MATLAB ++.

+1

All Articles