Chciałbym rozwiązać system liniowyAx = b w sposób liniowy najmniejszych kwadratów, uzyskując w ten sposóbx. MacierzeA, x ib zawierają elementy, które są liczbami złożonymi.

MatrycaA ma wymiaryn przezn, iA jest kwadratową matrycą, która jest również niższa trójkątna. Wektoryb ix mieć długościn. Istnieje tyle niewiadomych, ile jest równań w tym systemie, ale od tego czasub jest wektorem wypełnionym rzeczywistymi zmierzonymi „danymi”, podejrzewam, że byłoby lepiej zrobić to w sposób liniowy najmniejszych kwadratów.

Szukam algorytmu, który sprawnie rozwiąże ten system w sposób LLS, używając być może rzadkiej struktury danych macierzy dla macierzy o niższych trójkątachA.

Być może istnieje biblioteka C / C ++ z takim algorytmem już dostępnym? (Podejrzewam, że najlepiej jest użyć biblioteki ze względu na zoptymalizowany kod.) Rozglądając się w bibliotece macierzy Eigen, wydaje się, że dekompozycja SVD może być użyta do rozwiązania układu równań w sposób LLS (link do dokumentacji Eigen). Jak jednak pracować z liczbami złożonymi w Eigen?

Wydaje się, że biblioteka Eigen współpracuje z SVD, a następnie wykorzystuje ją do rozwiązywania LLS.

Oto fragment kodu pokazujący, co chciałbym zrobić:

#include <iostream>
#include <Eigen/Dense>
#include <complex>

using namespace Eigen;

int main()

{

    // I would like to assign complex numbers
    // to A and b

    /*
    MatrixXcd A(4, 4);
    A(0,0) = std::complex(3,5);     // Compiler error occurs here
    A(1,0) = std::complex(4,4);
    A(1,1) = std::complex(5,3);
    A(2,0) = std::complex(2,2);
    A(2,1) = std::complex(3,3);
    A(2,2) = std::complex(4,4);
    A(3,0) = std::complex(5,3);
    A(3,1) = std::complex(2,4);
    A(3,2) = std::complex(4,3);
    A(3,3) = std::complex(2,4);
    */

    // The following code is taken from:
    // http://eigen.tuxfamily.org/dox/TutorialLinearAlgebra.html#TutorialLinAlgLeastsquares

    // This is what I want to do, but with complex numbers
    // and with A as lower triangular

    MatrixXf A = MatrixXf::Random(3, 3);
    std::cout << "Here is the matrix A:\n" << A << std::endl;
    VectorXf b = VectorXf::Random(3);
    std::cout << "Here is the right hand side b:\n" << b << std::endl;
    std::cout << "The least-squares solution is:\n"
    << A.jacobiSvd(ComputeThinU | ComputeThinV).solve(b) << std::endl;
}// end

Oto błąd kompilatora:

 error: missing template arguments before '(' token

AKTUALIZACJA

Oto zaktualizowany program pokazujący, jak radzić sobie z rozwiązaniami LLS za pomocą Eigen. Ten kod rzeczywiście kompiluje się poprawnie.

#include <iostream>

#include <Eigen/Dense>

#include <complex>


using namespace Eigen;


int main()

{

    MatrixXcd A(4, 4);
    A(0,0) = std::complex<double>(3,5);
    A(1,0) = std::complex<double>(4,4);
    A(1,1) = std::complex<double>(5,3);
    A(2,0) = std::complex<double>(2,2);
    A(2,1) = std::complex<double>(3,3);
    A(2,2) = std::complex<double>(4,4);
    A(3,0) = std::complex<double>(5,3);
    A(3,1) = std::complex<double>(2,4);
    A(3,2) = std::complex<double>(4,3);
    A(3,3) = std::complex<double>(2,4);

    VectorXcd b(4);
    b(0) = std::complex<double>(3,5);
    b(1) = std::complex<double>(2,0);
    b(2) = std::complex<double>(8,2);
    b(3) = std::complex<double>(4,8);

        std::cout << "Here is the A matrix:" << std::endl;
    std::cout << A << std::endl;

        std::cout << "Here is the b vector:" << std::endl;
        std::cout << b << std::endl;

    std::cout << "The least-squares solution is:\n"

        << A.jacobiSvd(ComputeThinU | ComputeThinV).solve(b) << std::endl;


}// end