用线性最小二乘法求解带有复杂元素和下三角正方形A矩阵的系统Ax = b [英] Solving system Ax=b in linear least squares fashion with complex elements and lower-triangular square A matrix

问题描述

我想以线性最小二乘法求解线性系统Ax = b，从而获得x.矩阵A，x和b包含复数元素.

I would like to solve the linear system Ax = b in a linear least squares fashion, thereby obtaining x. Matrices A, x and b contain elements that are complex numbers.

矩阵A的尺寸为n乘以n，并且A是方形矩阵，也是较低的三角形.向量b和x的长度为n.这个系统中的未知数与方程式一样多，但是由于b是一个充满实际测得的数据"的向量，我怀疑最好以线性最小二乘法进行.

Matrix A has dimensions of n by n, and A is a square matrix that is also lower triangular. Vectors b and x have lengths of n. There are as many unknowns as there are equations in this system, but since b is a vector filled with actual measured "data", I suspect that it would be better to do this in a linear least squares fashion.

我正在寻找一种算法，该算法将以稀疏矩阵数据结构用于下三角矩阵A，从而以LLS方式有效地解决该系统.

I am looking for an algorithm that will efficiently solve this system in a LLS fashion, using perhaps a sparse matrix data structure for lower-triangular matrix A.

也许有一个C/C ++库已经可以使用这种算法了? (由于优化的代码，我怀疑最好使用一个库.)在本征矩阵库中四处查看，似乎可以使用SVD分解以LLS方式求解方程组(

Perhaps there is a C/C++ library with such an algorithm already available? (I suspect that it is best to use a library due to optimized code.) Looking around in the Eigen matrix library, it appears that SVD decomposition can be used to solve a system of equations in a LLS fashion (link to Eigen documentation). However, how do I work with complex numbers in Eigen?

Eigen库似乎与SVD一起使用，然后将其用于LLS解决.

It appears that the Eigen library works with the SVD, and then uses this for LLS solving.

这是一个代码片段，展示了我想做的事情:

Here is a code snippet demonstrating what I would like to do:

#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

这是编译器错误:

 error: missing template arguments before '(' token

---

更新

这里是更新的程序，显示了如何使用Eigen处理LLS求解.这段代码确实可以正确编译.

Here is an updated program showing how to deal with the LLS solving using Eigen. This code does indeed compile correctly.

#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

推荐答案

由于std::complex是模板类，因此使用std::complex(1,1);进行初始化时，编译器不知道它是什么类型.

Since std::complex is a template class, and you init with std::complex(1,1); the compiler doesn't know what type it is.

改为使用std::complex<double>(1, 1);.

这篇关于用线性最小二乘法求解带有复杂元素和下三角正方形A矩阵的系统Ax = b的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！