Как реализовать левое матричное деление на C ++, используя gsl - PullRequest
Новая
       22

Как реализовать левое матричное деление на C ++, используя gsl

5 голосов
/ 31 октября 2011

Я пытаюсь портировать программу 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) и не очень много знаю о математике, хотя быквадратная подгонка или матричная операция, может кто-нибудь сказать мне, как это сделать в 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;
 }

Ответы [ 2 ]

4 голосов
/ 31 октября 2011

В дополнение к тому, который вы дали, быстрый поиск выявил другие 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 , который возвращает минимальное решение по методу наименьших квадратов, в дополнение к ряду других итерационных методов для решения систем линейных уравнений.

1 голос
/ 31 октября 2011

Я не уверен, что понимаю ваш вопрос, но если вы уже нашли свое решение с помощью MATLAB, вы можете рассмотреть возможность использования MATLAB Coder , который автоматически переводит ваш код MATLAB в C ++.

...