【问题标题】:Solving system Ax=b in linear least squares fashion with complex elements and lower-triangular square A matrix以线性最小二乘方式求解系统 Ax=b,具有复杂元素和下三角方形 A 矩阵
【发布时间】:2012-11-24 16:46:27
【问题描述】:

我想以线性最小二乘方式求解线性系统Ax = b,从而得到x。矩阵A、x 和b 包含复数元素。

矩阵A 的尺寸为n 乘以n,而A 是一个也是下三角形的方阵。向量b 和x 的长度为n。这个系统中的未知数和方程一样多,但由于b 是一个充满实际测量“数据”的向量,我怀疑最好以线性最小二乘方式进行。

我正在寻找一种能够以 LLS 方式有效解决该系统的算法,可能使用稀疏矩阵数据结构来处理下三角矩阵 A。

也许已经有一个具有这种算法的 C/C++ 库? (由于代码优化,我怀疑最好使用库。)环顾 Eigen 矩阵库,似乎 SVD 分解可用于以 LLS 方式求解方程组 (link to Eigen documentation)。但是,如何在 Eigen 中处理复数?

Eigen 库似乎与 SVD 一起使用,然后将其用于 LLS 求解。


这是一个代码 sn-p 演示我想做的事情:

#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 求解。这段代码确实可以正确编译。

#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

【问题讨论】:

  • 那个错误指的是哪一行?
  • @RodyOldenhuis:行指A(0,0) = std::complex(3,5);,对象类型赋值std::complex。
  • 你能展示一下什么是 MatrixXcd 吗?
  • 等等,你的矩阵已经是下三角矩阵了?为什么不直接使用forward and back substitution?
  • 你不必使用 std::complex(2,4);代替?

标签: c++ linear-algebra eigen


【解决方案1】:

由于std::complex 是一个模板类,而您使用std::complex(1,1); 进行初始化,编译器不知道它是什么类型。

请改用std::complex&lt;double&gt;(1, 1);。

【讨论】:

  • 谢谢;你说得对,编译器不知道对象的类型。
猜你喜欢
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 2016-08-13
  • 2018-08-06
  • 2014-04-05
  • 2012-12-15
  • 1970-01-01
  • 1970-01-01
相关资源
最近更新 更多