5

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

A is an m-by-n matrix with m is not equal to n and B is a column vector with m components.

And I want the result X = A\B is the solution in the least squares sense to the under- or overdetermined system of equations AX = B. In other words, X minimizes norm(A*X - B), the length of the vector AX - B. That means I want it has the same result as the A\B in MATLAB.

I want to implement this feature in GSL-GNU (GNU Science Library) and I don't know too much about math, least square fitting or matrix operation, can somebody told me how to do this in GSL? Or if implement them in GSL is too complicate, can someone suggest me a good open source C/C++ library that provides the above matrix operation?

Okay, I finally figure out by my self after spend another 5 hours on it.. But still thanks for the suggestions to my question.

Assuming we have a 5 * 2 matrix 

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

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

The solution to the A \ b would 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;
 }
Amro
  • 123,847
  • 25
  • 243
  • 454
Daniel Gao
  • 207
  • 1
  • 3
  • 9
  • 1
    If you don't know linear algebra and numerical analysis, you shouldn't take it upon yourself to be implementing something like this, because you don't how to watch out for things like numerical instability. Just use GNU Octave, a free MATLAB replacement that has a good C++ API. – user57368 Oct 31 '11 at 02:28
  • 1
    Just so you know, matrix division isn't something that's typically well defined. You usually try to find a matrix inverse. – Mike Bailey Oct 31 '11 at 02:59
  • @DanielGao: maybe you should have posted your code as an answer – Amro Oct 31 '11 at 17:40
  • related question: [System of linear equations in C++?](http://stackoverflow.com/questions/2474432/system-of-linear-equations-in-c) – Amro Oct 31 '11 at 23:06
  • 2
    @mike why would you want to invert a matrix? Very expensive. Also impossible for non square as here. Division is exactly what is needed here. – David Heffernan Oct 31 '11 at 23:10
  • @Amro I actually try to post it as an answer, but the system does not allow me to do so since I do not have enough permission – Daniel Gao Nov 01 '11 at 04:15

2 Answers2

4

In addition to the one you gave, a quick search revealed other GSL examples, one using QR decomposition, the other LU decomposition.

There exist other numeric libraries capable of solving linear systems (a basic functionality in every linear algebra library). For one, Armadillo offers a nice and readable interface:

#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;
}

Another good one is the Eigen library:

#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;
}

Now, one thing to remember is that MLDIVIDE is a super-charged function and has multiple execution paths. If the coefficient matrix A has some special structure, then it is exploited to obtain faster or more accurate result (can choose from substitution algorithm, LU and QR factorization, ..)

MATLAB also has PINV which returns the minimal norm least-squares solution, in addition to a number of other iterative methods for solving systems of linear equations.

Amro
  • 123,847
  • 25
  • 243
  • 454
  • Thanks for the nice answer. Can Armadillo library solve linear equation AX=b, where A is a square sparse matrix with size 10000 by 10000? Does it apply the same algorithm in Matlab or Scilab, in which if A is square and nonsingular, x=A\b (backslash is defined for left matrix division) is equivalent to x=inv(A)*b (but the computations are much cheaper)? – Timothy May 06 '14 at 18:50
  • @Jeremy: According to the [docs](http://arma.sourceforge.net/docs.html#SpMat), Armadillo's support for sparse matrices is still preliminary, and some operations are still not fully optimized or implemented. As I [previously explained here](http://stackoverflow.com/a/18553768), `mldivide` in MATLAB (`x=A\b`) uses UMFPACK library to implement direct solvers for sparse linear systems. Just remember that in general you should never numerically solve a system using matrix inverse. In your case (a general square non-symmetric matrix) you could use the sparse LU decomposition. – Amro May 06 '14 at 20:50
  • Eigen library seems to have a robust support of [solving sparse systems](http://eigen.tuxfamily.org/dox-devel/group__TopicSparseSystems.html). There are other libraries out there you could use: http://math.stackexchange.com/q/638712/133 – Amro May 06 '14 at 20:52
1

I'm not sure I understand your question, but if you've already found your solution using MATLAB, you may want to consider using MATLAB Coder, which automatically translates your MATLAB code into C++.

zanbri
  • 5,958
  • 2
  • 31
  • 41