Как реализовать деление левой матрицы на C++ с помощью gsl

Я пытаюсь перенести программу MATLAB на C++. И я хочу реализовать деление левой матрицы между матрицей A и вектором-столбцом B.

A - матрица m-by-n, причем m не равно n, а B - вектор-столбец с m компонентами.

И я хочу, чтобы результат X = A\B был решением в смысле наименьших квадратов недо- или переопределенной системы уравнений AX = B. Другими словами, X минимизирует norm(A*X - B), длину вектора AX - B. Это означает, что я хочу, чтобы он имел тот же результат, что и A\B в MATLAB.

Я хочу реализовать эту функцию в 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;
 }
5
задан Amro 31 October 2011 в 23:06
поделиться